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CHAPTER  1 

INTRODUCTION 

Nickel-cadmium  storage  catteries  have  been  in  use  for 
more  than  fifty  years:  since  World  War  II,  r.ew  design  princi¬ 
ples  have  been  introduced  which  have  greatly  extended  their 
usefulness.  Many  theories  and  investigations  have  been  put 
forward  concerning  the  reaction  mechanisms  of  both  cadmium  and 
nickel  electrodes.  Because  the  activity  of  the  nickel  electrode 
has  been  more  fully  characterized  than  that  of  the  cadmium 
electrode,  we  will  focus  our  attention  on  the  behavior  of  the 
cadmium  electrode. 

There  are  a  number  of  methods  used  for  the  fabrication 

of  cadmium  electrodes.  These  methods  are  usually  classified 

into  three  categories:  (a)  pressed  powder  or  pasted  typev ^ , 

(4  5  6  7  S ) 

(b)  impregnated  nickel  sinter  typev  ,-3t  ’  '  '  arid  (c)  electro- 

(9  10  11  12 ) 

chemical  impregnation  71  ’  ’  The  last  method,  however,  is 

simple  in  operation  and  can,  in  theory,  achieve  high  loading  of 
active  material  in  a  single  impregnation  cycle  at  the  lowest  cost. 

In  the  electrochemical  impregnation  process,  electro- 
chemically  active  material  is  introduced  into  the  porous  nickel 
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plaqe  by  electrochemical  precipitation.  The  plaque,  which  serves 
as  a  cathode,  is  submerged  into  an  electrolysis  cell  containing 
an  aqueous  solution  of  cadmium  nitrate  at  a  temperature  near 
the  boiling  point  of  the  solution.  The  initial  pH  value  of  the 
impregnation  solution  is  adjusted  to  around  3  by  titration  with 
nitric  acid.  (As  the  pH  value  of  the  bath  increases,  i.e., 
becomes  more  basic,  cadmium  hydroxide  begins  to  precipitate  in 
regions  other  than  within  the  pores  of  the  plaque,  and  active 
material  is  lost  in  the  bath,  which  is  not  what  we  want.) 

The  electrochemical  reaction  at  the  cathode  during  impre¬ 
gnation  could  either  be  the  liberation  of  hydrogen, 

2H20  +  2e"  - ►  H2  +  20H"  ,  (1) 

or  the  reduction  of  nitrate  ion, 


(2) 


3oth  reactions  produce  hydroxy  ions  which  then  either  diffuse 
away  from  the  electrode  porous  surface  or  react  with  the  cadmium 
ions  present  in  the  pores  to  form  cadmium  hydroxide.  The  cadmium 


H20  +  NO^  +  2e“ - »  N02"  +20H 

•  •  •  •  • 

2H£0  +  NO^"  +  2e”  -  HNC>2  +30H 

t  9  •  •  • 

6H20  +  ;i03"  +  8e~ - -  NH  +90H 
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hydroxide  then  deposits  within  the  pores  of  the  porous  plaque. 

The  precipitation  reaction  is: 

Cd+2  +  20H“  - *  Cd(  OH)  2 .  (3) 

There  is  the  possibility  of  forming  cadmium  according  to  the 
reduction  reaction: 

Cd(0H)2  +  2e~  - *-  Cd  +  20H".  (4) 

This  reaction  is  presumably  prevented  from  taking  place  by 
hydroxy  ion  generated  on  the  electrode  surface  by  reaction  (1) 
or  (2),  so  that  there  is  no  chance  for  cadmium  ions  to  diffuse 
through  the  solution  onto  the  surface  and  be  reduced  to  cadmium 
metal . 


Following  impregnation,  the  electrode  is  cathodically 
polarized  in  a  25  $  potassium  hydroxide  solution  with  the  impre¬ 
gnated  plate  connected  as  the  cathode  and  subjected  to  low  direct 

2 

anodic  current  (0.25  to  0.75  amp/in  ).  The  solution  is  maintained 
at  about  100°C  for  a  period  of  from  15  to  45  minutes  to  convert 
cadmium  hydroxide  to  cadmium  in  the  pores  of  the  plaque  by 
reaction  (4)  or  to  remove  any  nitrate  ion  or  other  undesirable 
reduction  products.  The  active  material  of  the  electrode  is  pre¬ 
sumed  to  be  both  cadmium  hydroxide  and  cadmium. 
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In  this  work,  the  electrochemical  impregnation  process 
for  the  fabrication  of  the  cadmium  electrode  was  examined  both 
experimentally  and  theoretically. 

Experimentally,  a  set  of  current-density- vs.-time  transients 
at  constant  applied  potential  were  obtained  for  a  nickel  micro¬ 
electrode.  The  current  density  was  due  to  the  occurrence  of  the 
electrode  reaction  (2)  and  precipitation  reaction  (3)  presented 
earlier.  For  the  purpose  of  identifying  the  electrochemical 
reactions,  another  set  of  current-density -vs.-time  transients  were 
obtained  using  an  electrolyte  which  contained  no  cadmium  ion. 

In  this  case,  the  precipitation  reaction  (3)  no  longer  took 
place.  The  reactions  at  the  electrode  were  the  only  reactions 
taking  place. 

The  transient  current  responsed  in  an  electrolyte  which 
contained  cadmium  ions  differed  significantly  from  that  in  an 
electrolyte  which  did  not  contain  the  cadmium  ions.  This  diff¬ 
erence  is  presumed  to  be  due  to  the  precipitation  reaction. 

The  current  obtained  in  the  presence  of  cadmium  ions  in  the 
electrolyte  at  any  time  was  higher  than  that  obtained  at  the 
same  applied  potential  in  the  absence  of  cadmium  ions.  This 
result  is  explainable  by  the  fact  that  the  cadmium  ions  in  the 
solution  consumed  the  hydroxy  ions  which  were  produced  by  the 
heterogeneous  reaction  (2)  and,  as  a  result,  promoted  the  hetero¬ 
geneous  electrode  reactions. 
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The  experimental  current -density-  vs.-time  curves  obtained 
in  the  potential-step  experiment  were  used  to  identify  both 
the  heterogeneous  electrochemical  and  the  homogeneous  preci¬ 
pitation  reaction  rate  expressions.  The  heterogeneous  rate  was 
identified  first  from  the  experimental  results  in  which  there 
was  no  cadmium  ion  in  the  electrolyte.  The  precipitation  rate 
was  then  identified  from  the  experimental  results  in  which  the 
cadmium  ions  were  present  in  the  electrolyte. 

A  set  of  differential  equations  which  described  the  trans¬ 
port  process  in  the  cadmium-ion-free  solution  was  solved  in  terms 
of  two  unknown  constants.  The  unknown  constants,  related  to  the 
heterogeneous  electrode  reactions,  were  identified  by  matching 
the  theoretical  results  with  the  results  of  the  experiment.  The 
precipitation  reaction  rate  constant  was  then  obtained  by  a 
similar  procedure. 
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CHAPTER  2 


LITERATURE  REVIEW 

Impregnation  of  nickel  sinters  with  cadmium  hydroxide, 
followed  by  electrochemical  reduction  to  cadmium,  is  the  most 
economical  method  by  which  to  fabricate  cadmium  electrodes. 

Almost  all  of  the  research  works  concerning  this  subject  may 
be  divided  into  three  types i  (a)  manufacturing  techniques  de¬ 
velopment,  (b)  impregnation  process  study  arid  (c)  charge-discharge 
study . 


During  recent  years,  manufacturing  techniques,  especially 
those  involving  electrochemical  impregnation,  have  employed 
widely  varying  operating  conditions^ .  There  are 
several  parameters  that  affect  the  electrochemical  impregnation 
process  such  as: 

1)  Temperature:  Beauchamp  emphasized  the  fact  that  high 

impregnation  bath  temperature  keeps  the  size  of  cadmium 
hydroxide  crystals  small,  and  thereby  obtaining  a  more  active 
electrode . 

2)  pH  value:  In  Beauchamp’s  process  ,  the  pH  value  was 

closely  controlled  by  buffering  the  cadmium  nitrate  solution 

/  A  '  \ 

with  sodium  nitrite.  Pickett  °'  in  a  similar  process, 
maintained  the  pH  value  of  the  cadmium  nitrate  bath  between 
3-5  by  titrating  the  solution  with  proper  amount;  of  nitric 
acid . 
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3)  Current  density:  The  current  density,  based  on  apparent 

( 14) 

electrode  surface,  is  another  important  variable.  Beauchamp 
used  D.  C.  current  of  4-8  amps/in^  while  Bulan  ^-5)  usecj 

2 

alternating  current  of  1.2-1. 6  amps/in  and  claims  that  the 
electrode  impregnated  by  alterning  current  technique  retained 
more  of  its  capacity  after  cycling  as  compared  to  electrodes 
fabricated  by  other  methods.  Pickett presents  another 
technique  similar  to  the  above  and  claims  that  his  technique 
yields  higher  loading  cadmium  electrodes  and  achieves  about 
85%  utilization  of  the  active  material. 

The  electrochemistry  of  the  Cd/Cd(0H)2  electrode  in  con¬ 
centrated  alkaline  solution,  typically  KOH  as  in  the  Ni-Ca 
battery,  has  been  studied  extensively  both  theoretically  and 
experimentally  ^ ^ .  Bennion,  et  al.^^  in  their 

theoretical  development,  present  two  models,  the  solution- 
diffusion  model  and  film  model,  in  considering  the  failure 
mechanism  in  the  performance  of  secondary  batteries  using  metal/ 
metal  salt  couples.  They  conclude  that  one  mode  of  failure  of 
the  metal/metal  oxide  porous  electrode  is  the  blocking  of  the 

pore  or  the  complete  covering  of  surface  by  product  crystallities . 

(18) 

Falk'  obtained  X-ray  diffraction  patterns  from  electrodes 
submerged  in  electrolyte  during  charge  and  discharge  by  means 
of  a  special  test  cell.  He  claims  that,  during  charge,  Cd(CH)2 

is  successfully  transformed  into  Ca  metal  in  the  cadmium  electrode, 
but  even  after  a  strong  overcharge,  this  transformation  is  not 
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complete,  as  indicated  by  the  X-ray  pattern  at  this  stage. 

Furthermore,  Cd(OH),  is  the  end  product  during  discharge  and 

there  is  no  evidence  of  the  presence  of  a  CdO  film  which  is 

believed  to  be  responsible  for  the  passivation  of  the  electrode. 

( 19 ) 

Gross  and  Glocking  summarized  the  results  of  recent  research 
by  characterizing  the  negative  cadmium  electrode  in  a  nickel- 
cadmium  cell  and  evaluating  some  of  the  published  results  about 
which  researchers  still  disagree  among  themselves,  such  as  the 
mechanism  for  passivation,  which  is  an  important  phenomena  that 
affects  the  capacity  of  the  nickel-cadmium  cell. 

In  comparison  to  the  extensive  studies  on  the  electro¬ 
chemistry  of  the  Cd/Cd(OH)2  electrode  in  KOH  solution,  there 
is  little  information  available  concerning  the  electrochemical 
reactions  taking  place  during  the  electrochemical  impregnation 
process  which  uses  cadmium  nitrate  solution  as  electrolyte. 

The  only  research  works  known  to  the  author  are  the  electro- 

(21  22 ) 

cnemical  impregnation  studies  using  cycling  voltammetry  ’ 

2  2 
They  used  a  0.013  cm  nickel  microelectrode  and  a  0.064  cm 

cadmium  microelectrode.  Some  important  conclusions  of  their 

studies  are  summarized  below: 

(1)  There  was  a  large  reduction  wave  at  -O.63V  vs  SGE( Saturated 
Calomel  Electrode)  which  appeared  only  during  the  initial 
negative  scan.  This  is  believed  to  be  the  reduction  of  a 

solution  species  to  liberate  0H~  ion  in  the  presence  of 

+2 

Cd  to  precipitate  cadmium  hydroxide  on  the  electrode 


surface . 
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(2)  A  second  reduction  process  was  observed  at  -0.62V  vs  SGE 
during  the  reverse  positive  scan.  This  is  due  to  the  reductive 
formation  of  metallic  cadmium  from  cadmium  oxide  or  hydroxide 
that  was  deposited  on  the  electrode  surface  during  the  pre¬ 
vious  forward  scan. 

(3)  Some  processes  (or  one  of  the  processes)  during  the  first 
cycle  effectively  passivated  the  electrode  to  further  electro¬ 
activity.  The  formation  of  the  gray  material  at  -O.63V  vs 
SCE,  discussed  in  (1),  was  the  most  probable  reason  for  the 
passivation  of  the  electrode. 
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CHAPTER  3 

THEORETICAL  APPROACH 

3 . 1  General 

In  voltammetry  at  constant  voltage,  the  current  through 
the  electrolytic  cell  is  measured  as  a  function  of  time  while 
the  potential  of  the  polarized  electrode  is  held  at  a  constant 
value.  The  current  is  affected  by  both  the  transfer  process  and 
the  electrochemical  and/or  chemical  reactions. 

A  plane  electrode  is  immersed  in  a  solution  which  contains 
a  substance  0.  At  the  imposed  potential,  0  is  reduced  to  another 
substance  R  at  the  electrode  surface,  according  to  the  following 
electrochemical  reaction: 

0  +  ne"  - ►  R  (5) 

The  current-density- vs.-time  curve  obtained  at  a  constant 
voltage  depends  on  both  the  kinetics  of  the  electrochemical 
reaction  and  the  rate  of  mass  transfer. 

Under  our  experimental  conditions,  diffusion  is  the  only 
mode  of  mass  transfer.  The  other  two  modes  of  mass  transfer, 
namely,  migration  and  convection,  are  insignificant.  Migration 

I 
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is  eliminated  by  the  presence  of  a  large  excess  of  supporting 
electrolyte  which  reduces  the  potential  field.  Convection  is 
avoided  due  to  the  conditions  that  (1)  the  solution  is  not 
stirred  and  (2)  the  duration  of  electrolysis  is  short,  so  that 
the  density  change  has  not  yet  become  a  significant  factor  and 
natural  convection  has  not  taken  place. 


The  dimension  of  the  cell,  as  compared  to  the  microelectrode 
used,  is  so  large  that  the  solution  may  be  regarded  as  extending 
to  infinity;  i.e.,  the  diffusion  process  is  semi-infinite. 
According  to  Fick's  first  law,  the  flux  of  the  substance  0  at 
a  distance  x  from  the  electrode  and  in  the  direction  perpen¬ 
dicular  to  the  electrode  is: 

^C0(x,t) 

N0  (x.t)  =  -  D0-“  (6) 

OX  • 

Conservation  of  the  species  0  leads  to  the  following  differential 
equation : 

?>Cn(xtt)  =  D  ^2c0(x,t) 

t  0  7)  x2  (?) 


where  Cq  is  the  concentration  of  substance  0,  which  is  a  function 
of  x  and  t.  Dq  is  the  diffusion  coefficient  of  substance  0, 
and  is  assumed  to  be  a  constant. 
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When  the  substance  0  is  reduced  at  a  very  high  rate,  by 
imposing  a  sufficiently  large  potential  step,  such  that  the 
concentration  of  this  substance  at  the  electrode  surface  is 
reduced  to  zero  soon  after  the  start  of  electrolysis ,  the 
boundary  condition  at  the  electrode/electrolyte  interface  for 
equation  (?)  is  CQ(0,t)=0  for  t>0.  Initially  before  the  electro¬ 
lysis,  the  solution  is  homogeneous,  and  the  concentration  Cq 
at  t=0  is  uniform  with  the  bulk  value.  The  other  boundary 
condition  is  that  CQ  approaches  its  bulk  value  when  x  is 
sufficiently  large. 

The  electrolysis  current  is  equal  to  the  product  of  the 
charge  involved  in  the  reduction  of  one  mole  of  substance  0  by 
the  flux  of  this  substance  at  the  electrode  surface.  Thus  the 
current  density  (current  per  unit  electrode  surface  area)  is: 

1=  nFN  0 ( 0 , t ) ,  (8) 

where  n  is  the  number  of  electrons  involved  in  the  reduction 
reaction  .of  substance  0,  F  is  the  Faraday  constant,  and  NQ(0,t) 
is  the  flux  of  substance  0  for  x=0. 

The  solution  of  the  above  problem  is^-^  : 

1=  nFD  *  cj  - x- ~  ” -  (9) 

U  U  t  2 

f 

where  Cq  is  the  bulk  value  of  substance  0. 
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The  current  density  represented  by  equation  (9)  is  called 
the  limiting  current  density  of  electrochemical  reaction  (5), 

This  is  the  maximum  current  density  which  can  be  achieved  by 
this  reaction.  It  was  found  that  if  nitrate  ion  is  the  only 
reducible  substance,  the  current  density  predicted  by  equation 
(9)  is  considerably  lower  than  the  experimental  value.  This 
contradition  is  resolved  if  one  realizes  that  under  the  experi¬ 
mental  conditions  water  is  also  reduced. 

At  any  other  less-negative  potential  the  current  will 
depend  on  the  electrode  kinetics.  Now  the  condition  thax  the 
reducible  substance  concentration  becomes  zero  after  the  appli¬ 
cation  of  the  potential  step  has  to  be  replaced  by  a  kinetics 
equation,  which  is  a  function  of  potential  and  the  concentrations 
of  surface  reactant  and  product.  No  general  solution  is  available 

except  for  the  case  when  the  reaction  rate  is  linearly  dependent 

,  .  ( 21 ) 
on  tne  concentrations  of  the  surface  reactant  and  product  J  . 

3.2  Description  of  the  Model 

The  real  porous  electrode  is  formed  by  sintering  nickel 
powder  on  a  nickel  grid  to  form  a  plaque.  The  plaque  serves  as 
the  cathode,  which  is  positioned  in  a  cadmium  nitrate  solution 
with  a  pH  value  between  3  to  5*  The  passage  of  the  current 
produces  hydroxy  ions  which  cause  the  cadmium  to  be  deposited 
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into  the  porous  pores.  The  transport  phenomena  associated 
with  electrochemical  reactions  in  a  porous  medium  is  very 
complicated.  The  approach  here  is  to  study  the  intrinsic 
reactions  which  take  place  inside  the  pores  in  a  simple 
flat  electrode.  The  plate  electrode  is  submerged  in  a 
semi-infinite  pool  of  electrolyte  so  that  the  mass  transfer 
problem  can  be  treated  as  a  one-dimensional  problem  in 
the  x-direction  (perpendicular  to  the  electrode  surface). 

In  the  solution  side,  there  is  a  thin  double-charge  layer 
near  the  electrode.  The  thickness  of  the  double  layer  is 

O 

on  the  order  of  10  to  100  A.  It  is  too  thin  to  be  probed 
adequately,  and  the  theory  of  the  diffuse  layer  is  a 
microscopic,  rather  than  a  macroscopic,  phenomena.  The 
mass  transfer  region  to  be  considered  here  is  thus  located 
outside  of  this  double  layer. 

The  electrodeposition  process  operates  on  the  principle 
that  the  hydroxy  ion,  which  is  generated  by  the  reduction 
process,  coprecipitates  with  the  cadmium  ion  in  the  solution 
to  form  cadmium  hydroxide  crystalloid.  The  approach  taken 
here  to  unravel  the  complicated  sequence  of  reactions  is 
to  form  a  model  which  describes  the  transport  processes 
and  reactions.  The  model  will  first  be  used  to  analyze 
the  experimental  results  obtained  under  the  condition  of 
no  pricipitation  reaction,  by  using  an  electrolytic  solution 
which  is  free  of  cadmium  ions.  The  heterogeneous  reduction 
reaction  will  be  determined.  The  transport  and  reacxion 
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problems  will  then  be  examined  in  the  case  when  the  solution 
contains  cadmium  ions.  The  rate  of  the  precipitation 
reaction  will  then  be  determined  by  matching  the  model's 
prediction  and  the  experimental  results. 


3 . 3  Formulation  of  Transport  Process  Equations 


Three  ionic  species  to  be  considered  in  the  electro- 

_  _  +2  . 

chemical  system  are  NO^ ,  GH  and  Cd  ions.  In  the  absence 
of  migration  and  convection,  the  flux  of  each  species  is 
expressed  as^^^s 


Ni  -  "  Di 


,  1  =  1, <2, 3  , 


where 


is  the  flux  of  species  i  (moles/cm  -sec ) , 
is  the  diffusion  coefficient  of  species  i  (cm  /sec), 
Ch  is  the  concentration  of  species  i  (moles/cur  ) , 

which  is  a  function  of  both  x  and  t,  and  i=l,2,3 

-  -  +2 

corresponding  to  NO^,  GK  and  Cd  ions. 

(24) 

The  equation  of  continuity  of  each  species  is 


?Ci  = 

"St 


d  Ni 


i=l,2,3  , 


where  rn  is  the  rate  of  production  of  species  i 
( raoles/cm  -sec )  by  precipitation  reaction,  which  occurs 
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only  in  the  presence  of  cadmium  ions. 


Substitution  of  equation  (10)  into  (11)  give: 


?>t 


2 

3  c. 


=  Di' 


2):. 


2 


i=l,2,3 


(12) 


It  is  assumed  that  the  precipitation  reaction  (3) 
can  be  treated  as  a  homogeneous  reaction  and  is  linearly 

proportional  to  the  degree  of  supersaturation.  Therefore, 

the  consumption  rate  of  cadmium  ion  due  to  this  reaction 

is  ; 


R 


3 


(13) 


where  k  is  the  homogeneous  reaction  rate  constant  (sec-^), 
K sp  is  the  solubility  constant  of  Cd(0H)2 

and  has  the  value  of  2xl0-^  (moles^/cm^)  ^5) . 

3y  the  stoichiometric  relation  of  equation  (3)>  the  con¬ 
sumption  rate  of  hydroxy  ion  due  to  the  same  reaction  is  : 


«2  =  2  R3 


-2k 


(c2)‘ 


(14) 


Substitution  of  equations  (13)  ana  (14)  into  equation  (12' 
leads  to  the  following  two  differential  equations  s 


7)t 


->c. 


-d° 2  =  2  k 


Ksp 


(c2)‘ 


(  If) 


2  { 

T'3 

D-. - k  C_ 

3}x2  I  3 


t  p  \  2 
v  2; 


For  NO"  ion, 


R1  =  o. 


Equation  (12)  for  NO^  ion  may  be  reduced  to 


2>C1  2Tci 

- —  =  n  - — 

■at  i  \  2 

O  X 


Only  two  of  the  three  differential  equations  need 

to  be  solved  for  the  concentration  profiles,  due  to  the 

(24) 

electroneutrality  condition  s 


2C3-C1-C2=0 


Substitution  of  equation  (18)  into  equation  (15) 


gives  s 


~2>  C2 
it 


Equations  (1?)  and  (19)  are  solved  simultaneously 

to  obtain  the  concentration  profiles  of  ar.d  C~,  i.e., 

-  .  +2 

i.'O^  and  OH  ions.  The  profile  can  then  be  obtained  by 
equation  (18)  from  the  known  profile^  of  and  CU . 


A  set  of  differential  equations,  which  are  simpler 
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than  equations  (17)  and  (19).  ars  obtained  when  the  elec 
trolyte  contains  no  cadmium  ion.  In  the  absence  of  cad¬ 
mium  ion,  the  precipitation  reaction  (3)  can  not  occur; 
thus,  FL  terms  are  all  equal  to  zero.  Equation  (12)  is 
then  reduced  to i 


7>Ci 

at 


i=l i 2  , 


(20) 


in  which  i=l,2  corresponds 


to  N0^ 


and  OH 


ions . 


Equation  (20)  is  written  for  and  C^s 


DC1 

7>t 


= 

l>t 


(21) 


(22) 


The  other  condition,  the  concentrations  of  various 
ions  at  the  electrode/electrolyte  interface  after  the 
start  of  the  current  flow,  has  to  be  described  by  the 
electrode  kinetics  of  the  reduction  reaction  of  N0“  ions 


3.^  Electrode  Kinetics 


The  possible  reduction  reactions  are  numerous,  as 
shown  in  equation  (2).  Since  in  e*..ch  reaction,  OH"  ion 
is  generated,  we  will  represent  the  production  of  OH" 
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ions  as  the  following  electrochemical  reaction : 

2H20  +  N 0 ^  T  2e“ - *  HN02  +  30H".  (23) 

Because  there  are  two  electrons  involved  in  the  reaction, 
it  is  likely  not  an  elementary  step.  Let's  assume  the 
reaction  comprises  two  elementary  steps,  each  step  involving 
the  transfer  of  one  electron,  and  a  dissociation  reaction: 


(slow) ,  (24a) 

(fast),  (24b) 

(fast),  (25) 


where  and  Kc2  are  cathodic  reaction  constants,  and  K  ^ 
and  Ka2  are  anodic  reaction  constants,  respectively.  is 
the  dissociation  reaction  constant.  (OH) 2  ion  is  an  unstable 
intermediate  ion  which  immediately  enters  the  next  reaction  (24b/ 
The  rates  of  these  elementary  steps  should  always  be  propor¬ 
tional  to  each  other.  For  the  elementary  steps  listed  above, 
reaction  (24b)  should  occur  once  every  time  reaction  (22a)  occurs 


Let  I.  and  I,  denote  the  cathodic  current  densities 
of  reactions  (24a)  and  (24b),  and  r,  and  r2  denote  the 
reaction  rates  of  reactions  (24a)  and  (24'c),  respectively. 
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Then,  tne  relationship  between  r.  ,  I-  and  surface  ion 


concentrations  of  N cZ,  (CK)2  ,  which  app 

( 2  ^ ) 

step  (24a),  may  be  written  as 


ear  in  e-emer.tar; 


-  _Ll 

'l"  nF 


(n=l ; 


3  Kal  (CNC;](C(0H)-)  fJ 

W  W 


( :-/ 3/ :? 


..  r„o  I  fro  1  o„^  r  f 

cl  °ho:  CH90  e'‘?  "  " 

u  JJ  u  ^  L  *" 


r  a?  „ 


Similarly,  for  reaction  (24b),  one  has: 


(n=l ) 


^K-r«pp^v] 

Ko2  [C(0K)-]  ex?[-^Lvj' 


where  (3,  and  ^  are  symmetry  factors  representing  the  fractions 
of  applied  potential  V,  which  promotes  the  cathodic  reaction. 
C.°.~-  ,  and  Z%,-  are  the  surface  ion  concentrations  of 

4'0y  ( OH ) ^  and  OH”  on  the  electrode  surface. 


Since  reactions  (24a)  and  (2^b)  occur  at  the  same  rate, 
we  have : 


1=  1^  f  I,  s  2I2  =  21.  , 


2 


21 


where  I  is  overall  current  density.  It  is  further  assumed 
that  reaction  (24b)  is  fast,  and  is  essentially  in  equilibrium. 
From  equation  (2?), 


-r 

xc2  exp  ['-fr7] 


and  = 
r 


=r  =  overall  reaction  rate. 


Substituting  equation  (28)  into  equation  (25)  and 
combining  equation  (29)  gives: 


_  KalKa2  (  ro 


[C°°-][C0OH-]2  exp£ 


(2-A)? 


xd(ckJfS2o]-p[-^v]  . 


From  reaction  (25)  one  may  write 


Kd  = 


(GHN02][g0H~] 

[cno“][ch2oJ 


Substitution  cf  the  above  expression  into  equation  ijC)  leads 
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ions  and  electrolyte  containing  no  cadmium  ions,  two 
and  four  boundary  conditions  are  needed  to  solve  the 
equations  (17)  and  (19)  or  (21)  and  (22). 


initial 
differ en 


At  t=0,  the  concentration  of  each  ion  is  equal  to  its 
bulk  value  at  each  point  along  the  x-direction.  The  initial 
conditions  are  as  follows: 

1#  at  t=0,  C^=  for  all  x, 

2.  at  t=0 ,  C2=  C 2  for  all  x, 

where  the  superscript  "c"  denotes  the  bulk  condition. 


The  first  boundary  condition  is  the  consequence  of 
the  fact  that  the  flux  at  the  interface  has  to  be  equal  to 
the  rate  of  heterogeneous  electrode  reaction,  i.e.: 

xm0  *  K2Cc°)3  -  XlC°?l  <32) 

where  superscript  "o"  denotes  the  surface  condition.  To 
satisfy  the  stoichiometric  relation  of  equation  (23),  one 
has,  at  x=0 : 


*  ^ 

,  >:1 

-  - 

f  \ 

>C2 

dx 

1 

o 

II 

X 

^2  -v 

3x 

^  J 

k  J 

This  is  the  second  boundary  condition.  The  ct.-.er  two 
boundary  conditions  are: 
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3  •  at  x=oo  ,  C^=  C° , 

•u 

4.  at  x=co  ,  C2=  C2. 

That  is,  the  concentrations  far  away  from  the  electrode/ 
electrolyte  interface  do  not  change. 


The  flux  in  equation  (32)  is  proportional  to  the  curren 
density  as  a  function  of  time: 
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CHAPTER  lr 

EXPERIMENTAL  WORKS 

A  potential-step  was  generated  from  a  function 
generator  and  was  applied  to  the  working  electrode/elec- 
trolyte  interface  in  the  electrochemical  cell.  The  ca¬ 
thodic  current  was  thus  sampled  and  stored  in  memory  at 
a  constant  time  interval  by  a  microcomputer.  The  resulting 
current/time  transient  was  then  plotted  by  a  x-y  recor¬ 
der  or  printed  out  by  a  line  printer  by  recalling  these 
current  data  from  the  memory  of  the  microcomputer. 

4. 1  The  Electrochemical  Cell 


The  current  transient  under  a  potential  step  con¬ 
dition  was  obtained  in  a  three-electrode  cell.  The  cell 
and  electrodes  configuration  is  sho-wn  in  Figure  1.  The 
working  electrode  was  a  nickel  microelectrode.  It  was  made 

of  a  nickel  wire  which  was  pressure-fitted  into  a  teflon 

.  2 
cylinder.  The  flat,  exposed  surface  was  0.013  cm4" .  A 

saturated  calomel  electrode  (SCE)  and  a  cadmium  bar  were 
used  as  the  reference  and  che  counter  electrodes,  res¬ 
pectively  . 


The  cell  was  filled  with  electrolyte  of  either 
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cadmium  nitrate  solution  or  potassium  nitrate  solution. 

The  pK  value  of  the  solution  was  controlled.  Potassium 
chloride  was  used  as  the  supporting  electrolyte.  The  electro¬ 
lytic  solution  was  not  stirred  during  the  experiment,  so  that 
the  conditions  of  semi-infinite  linear  diffusion  were  maintained. 

Two  experimental  parameters  were  varied: 

1)  the  ion  concentration  of  the  solutions  and 

2)  the  magnitude  of  the  potential-step  applied  to  the 
interface. 

4.2  The  Instruments  and  Electrical  Setup 

The- block  diagram  for  the  experiment  setup  is  shewn  in 
Figure  2.  The  instruments  shown  in  Figure  2  contain  the 
following  equipment: 

a)  A  Princeton  Applied  Research  (PAR)  Model  1'~’3  Potentiostat/ 
Galvan os tat :  In  the  experiment,  the  operating  mode  was  set 
at  CONTROL  E. ,  which  means  that  the  instrument  functioned 
as  a  potentiostat.  The  current  was  measured  while  keeping 
the  potential  of  the  working  electrode  constant.  This  was 
accomplished  by  setting  the  counter  electrode  to  the 
required  level  so  that  the  'working  electrode  potential  was 
at  a  programmed  potential  with  respect  to  the  reference 
electrode.  The  instrument  is  provided  with  a  Me  del  l"o 
Current- to -hoi tags  Converter  which  is  the  "basic"  piug-in 
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module.  The  M  176  provides  a  do  voltage  output  ’./hi oh  is 
proportional  to  the  current.  The  cell  current  can  be  eitner 
displayed- on  the  panel,  by  a  x-y  recorder,  or  monitored  by 
the  microcomputer.  Connection  to  the  cell  was  made  through 
the  external  cable »  The  counter  electrode  was  connected  to 
the  red  clip,  the  working  electrode  to  the  green  clip  and 
the  reference  electrode  to  the  Electrometer  Probe  which  has 
a  very  high  impedance,  thereby  insuring  us  that  the  current 
was  flowing  to  the  reference  electrode, 

b)  A  PAR  Model  175  Function  Generator!  This  is  a  programmable 
waveform  generator  which  has  two  operating  modes,  SWEEP  and 
PULSE.  The  former  generates' a  sequence  of  triangular  waves, 
while  the  latter  generates  a  sequence  of  step  function.  Eor 
the  potential  step  experiment,  the  operating  mode  was  set- a 
PULSE  mode  and  the  pulse  width  selector  at  STEP  position. 
The  magnitude  of  the  step- function  was  set  by  the  setting  o 
the  3  potential  (upper  limit)  in  the  POTENTIAL  section  on 
the  front  panel  of  this  instrument. 

c)  An  Intel  CPU  S080A  Bases  Microcomputer  ( cU-K  ?,A The  micro 
computer  was  the  central  part  of  the  experimental  setup. 

The  functions  of  the  microcomputer  in  the  experiment  were: 

(1)  to  trigger  the  Function-' Generator  which  activated  the 
potential  step  presented  earlier, 

(2)  ...to  sample  the  current  outputs  from  the  M  17c  Converter 

-  -and  store -them  in- the- memory .  — 

(3)  to  convert  these  digital  data  back  into  analog  signals 
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and  then  plot  them  on  the  x-y  recorder. 

(4)  to  print  the  digital  data  on  the  line  printer. 

(5)  to  store  these  digital  data  in  the  tapes. 

The  peripheral  equipment  of  the  microcomputer  include  a 
CRT  (TV  screen),  the  keyboard,  a  cassette  recorder  and  a 
line  printer. 

d)  A  Houston  Model  RE  0074  x-y  Recorder:  The  Recorder  receives 
the  analog  outputs  from  the  DAC  (Digital-Analog  Converter) 
and  plots  them  on  the  graph  paper. 

4.3  The  Computer  Program 

A  computer  program  ’written  in  BASIC  language  -was  used  no 
carry  out  the  experiment.  The  program  consists  of  a  main  program 
and  a  machine-language  subroutine  called  "MUG".  This  subroutine 
includes  an  execution  statement  which  generates  a  trigger  signal 
to  activate  the  Function  Generator.  A  potential  step  was  then 
immediately  applied  to  the  working  electrode,  and  thereafter 
the  computer  started  to  sample  and  store  the  current  output  at 
a  fixed  time  interval.  The  analog  signal  was  converted  into  a 
digital  datum  by  a  12  bits  ADC  (Analog-Digital  Converter)  and 
then  stored  in  the  memory.  When  the  number  of  samples  reached 
a  preset  value,  the  sampling  routine  was  terminated.  The  digital 
sigr.als  stored  in  the  memory  were  binary  numbers.  These  data 
were  converted  into  decimal  numbers  and  then  converted  into 
analog  signals,,  and  finally  plotted  on  the  x-y  recorder  and 
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printed  on  the  line  printer.  This  was  executed  by  the  main 
program.  The  main  program  and  the  "MUG"  subroutine  are 
presented  fully  in  Appendix  A1  and  A2 . 

4.4  Experimental  Conditions  and  Procedures 


In  the  experiment,  the  initial  pH  value  of  the  electro¬ 
lyte  was  maintained  at  3*G,  which  was  the  condition  used  in  the 
actual  electrochemical  impregnation  process.  This  was 
accomplished  by  titrating  the  electrolytic  solution  with  a 
concentrated  KC1  solution.  The  concentrations  of  nitrate 
ion  used  were  0.005M,  0.002574  and  0.00125M  prepared  from  either 
the  KNC^  or  Cd(NO^)  solution.  The  cadmium  ion  concentrations 
were  0.0025M,  0.00125M  and  0.00Q625M,  All  the  electrolytic 
solutions  contained  0.5M  KCi  as  the  supporting  electrolyte. 

It  was  .found  later  that  the  electrolysis  of  water  made 
significant  contribution  to  the  total  current,  because  the 
current  obtained  was  significantly  higher  than  the  limiting 
current.  In  order  to  obtain  the  current  due  only  to  the  reduction 
of  nitrate  ions,  several  .runs  with  solutions  which  contained  no  . 
nitrate  ion  were  conducted.  The  solutions  used  for  this  purpose 
were  either  water  or  C.0G25M  CdCl^  solution.  Beth  solutions 
contained  0.5 M  KCi  as  the  supporting  electrolyte.  This  current 
was  subtracted  from  the  current  obtained  in  t he  ores er.ee  of 
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nitrate  ions  in  the  solution.  The  resulting  current  is  the 
faradic  current  of  the  reduction  of  nitrate  ions. 

The  magnitudes  of  the  potential-step  used  in  this 
study  were  -0740V,  -0 . 60V  and  -0 . 80V  versus  the  equilibrium 
potential  which  was  at  approximately  -G.35^V  vs  SCE, 

The  experimental  procedures  are  listed  blows 

(1)  The  electrical  circuit  was  set  up  as  shown  in  Figure  2, 
except  that  the  cell  was  disconnected  from  the  Potent iostat/ 
Galvanostat  by  setting  the  cell  selector  on  the  instrument 
to  the  OFF  position.  A  proper  current  range  was  then 
chosen  from  the  Current  Range  Switch  to  make  sure  that 

the  current  output  would  not  be  overloaded.  The  Function 
Generator  was  initialized  by  depressing  the  INITIAL 
control  pushbuttom.  The  initial  potential  was  set  at 
zero  volts. 

(2)  A  known  volume  of  electrolytic  solution  and  an  equal 
volume  of  supporting  electrolyte,  KC1,  were  added  to  the 
cell  and  titrated  with  HC1  solution  to  a  pH  value  of  3*0. 

The  solution  was  then  deaerated  by  bubbling  purified 
nitrogen  through  the  stirred  solution  for  about  ten 
minutes.  The  gas  continued  to  pass  above  the  solution 
during  one  experiment. 

(3)  After  the  working  electrode  surface  was  polisr.sd  by  a 

3/0  alumina  paper,  it  was  rinsed  thoroughly  with  distilled 
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and  deionized  water,  and  then  was  positioned  in  the  cell. 

The  electrodes  (the  working  electrode,  counter  electrode 
and  reference  electrode)  were  then  connected  to  the 
Potentiostat/Galvanostat  as  described  in  section  4. 3 • 

After  the  electrodes  were  correctly  connected,  the  cell 
selector  was  switched  to  the  EXTERNAL  CELL  position.  The 
rest  potential  of  the  solution  could  be  determined  by 

turning  the  knob  in  the  APPLIED  POTENTIAL/ CURRENT  section 
of  the  front-panel  of  the  Potentiostat/Galvanostat  so 
that  the  current  reading  displyed  on  the  front  panel 
meter  was  zero.  The  applied  potential  was  then  the  rest 
potential. 

(4)  The  upper  limit  (B  potential)  on  the  Function  Generator 
was  set.  The  value  of  the  -3  potential  plus  the  rest 
potential  was  the  total  potential  applied  to  the  electror 
chemical  cell. 

(5)  The  experiment  was  initiated  by  running  the  computer 
program.  The'  current/time  data  were  stored  in  the  computer’s 
memory  and  -then  the  results  were  converted  to  line 

printer  output. 

4. 5  Experimental  Results 

Two  sets,  a  total  of  twenty-two  runs,  were  made.  The 
first  set  of'  runs  was  conducted  in  solutions  which  contained 
no  cadmium  ion.  The  second  set  of  runs  was  made  in  solutions 
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which  contained  cadmium  ions.  For  each  set  of  runs,  clank  runs 
were  made  in  the  solutions  without  the  presence  of  nitrate 
ion  for  the  purpose  of  determining  the  background  current 
aue'to  the  water  decomposition  reaction.  The  solution  for 
blank  runs  for-the  first  set  was. water,  while  the  solution 
for  blank  runs  for  the  second  set  was  made  by  using  the 
0.0025M  CdC^  solution. 


The  current-density-vs . -time  data  printed  by  the  H  14 
Line  Printer  are  shown  in  Table  3.1  through  Table  B.22  in 
appendix  3.  The  measured  current  was  converted  to  the  current 
density  which  appears  in  those  tables.  This  is  calculated 
by  dividing  the  current  by  the  electrode  surface  area 
0.013cm2. 


Tables  3.1  through  3.9  are  the  first  data  set  which  used 
potassium  nitrate  as  the  electrolytic  solution,  while  Tables 
3.10  through  3,12  are  the  background  currents  of  water  used 
to  subtract  from  the  first  data  set. 


Tables  3.13  through  3.21  are  the  second  data  set  which 
used  cadmium  nitrate  as  the  electrolytic  solution,  while 
Table  3.22'  is  the  background  current  of  water  in  tr.e  presence 
of  cadmium  ions. 


A  series  of  current-density-vs . -time  curves  with  various 
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electrolytic  solution  concentrations  and  applied  overpotentials 
are  shown  in  Figures  3  through  6.  In  these  figures,  the 
current  due  to  the  electrolysis  of  water  has  been  elirainited 
by  using  the  routine  tabulated  in  Table  1. 

Figures  3  through  5  are  the  current-density-vs . -time 
curves  obtained  in  the  solutions  which  did  not  contain  cadmium 
ion.  Figure  3  shows  the  current-density-vs . -time  curves 
obtained  in  a  solution  containing  0.Q05M  nitrate  ion  and  at 
potentials  -0.40V,  -0.60 V  and  -0.80V  from  the  equilibrium 
potential.  Figure  4  shows  the  current-der.sity-vs , -time  curves 
obtained  in  a  solution  containing  0.0025M  nitrate  ion,  and  at 
potentials  -0.40V,  -0.60V  and  -0.807  from  the  equilibrium 
potential.  Figure  5  shows  the  current-density-vs . -time  curves 
obtained  in  a  solution  containing  0.00125M  nitrate  ion  ax  -0.40V, 
-0.60V  and  -0„80V_from  the  equilibrium  potential.  Due  to  the 
lack  of  blank  runs  with  O.OOI25M  and  0.000625M  cadmium  ion  con¬ 
centrations  and  -0.60V  and  -0.80V  from  the  equilibrium  potential, 
only  one  current-density-vs . -time  curve  was  obtained  with  0.0G5M 
nitrate  ion  and  0.0025M  cadmium  ion  concentration  at  the  poten¬ 
tial  of  -0.40V  from  the  equilibrium  potential.  This  is  shown 
in  Figure  6. 

4. 6  Discussion  of  Experimental  Results 


Some  characteristics  can  be  seen  from  the  curves  shown 
in  Figures  3  through  6. 
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Table  Is  Evaluation  of  Currer.t-Density-vs .  -Time  lata 
Plotted  in  Figures  3  Through  6. 


— 

Figure  is 
plotted 

by  subtracting  the 
corresponding  current 
density  at  each  time 

uoint  in 

from  that 

of  in 

(1) 

. 

Table  3.1 

Table  3.10 

Figure  3 

(2) 

Table  3.2 

Table  3.' 11 

(3) 

Table  3.3 

Table  3. 12 

(1) 

Table  3.4 

Table  3.10 

Figure  4 

(2) 

Table  3.5 

Table  3.11 

(3) 

Table  3.6 

Table  B.12 

i 


(1) 


Table  3.7 


Table  3.10 


Firstly,  the  current  density  decays  with  tine  in  a  ] 

hyperbolic  fashion.  This  is  the  result  of  depletion  of  the  j 

iN'O^  ions  near  the  electrode  surface.  Theoretically,  there  j 

is  an  initial  sharp  rise  of  current  associated  with  double-layer 
(27) 

charging  ,  due  to  the  application  of  a  potential-step .  This  was  \ 

not  observed  in  the  experiment,  however,  because  the  tine  for  j 

charging  is  usually  very  short  as  compared  to  the  electrolysis  1 

time. 

;i 

i! 

Secondly,  the  sharp  current  drops  in  the  beginning  indicate  j 

that  very  sharp  concentration  gradients  are  established  for  the  i 

:  I 

nitrate  ions  and  the  hydroxy  ions  as  soon  as  the  surface 

reaction  (23)  takes  place.  The  concentration  gradients  are  ! 

i 

gradually  decreased  due  to  the  diffusion  layer  extending  in  ^ 

the  direction  of  decreasing  concentration  gradients.  The  i 

result  is  that  the  decay  of  the  current  density  is  not  as  fast  j 

as  in  the  beginning. 


The  curves  for  the  case  of  potential  at  -0.80V  from  the 
equilibrium  potential  and  solutions  containing  0.00251(1  and 
0.00125M  nitrate  ion,  as  shown  in  the  upper-most  curves  in  Figures  - 
and  5»  respectively,  are  somewhat  different  from  tne  others. 

These  two  curves  do  not  follow  the  same  trends  as  their  lower 
potential  counterparts.  It  is  suggested  that  t he  reactant, 
nitrate  ion,  near  the  electrode  surface,  is  used  up  in  a  very 
short  time  and  another  reaction  must  be  taxing  place.  This 


behavior  is  the  mos*  pronounced  as  shown  in  the  upper-most 
curve  in  Figure  5*  A  maximum  current  is  observed  at  time 
equal  to  about  0.5  msec.  This  curve  thus  strongly  supports 
our  explanation  that  some  other  reaction  is  taking  place  at 
that  time.  On  the  other  hand,  the  absence  of  this  behavior 
in  the  curve  for  the  case  of  0.005M  nitrate  ion  shown  in  the 
upper-most  curve  in  Figure  3  suggests  that  the  nixrate  ion 
concentration  near  the  electrode  surface  is  not  zero. 

Two  parameters  were  varied  in  the  experiment,  i.e., 
the  concentration  of  the  electrolyte  solutions  and  the  value 
of  the  potential  from  the  equilibrium  potential  (rest  potential). 
The  changing  of  either  the  bulk  concentration  of  the  electrolyte 
solutions  or  the  overpotential  that  is  applied  to  the  working 
electrode  will  change  the  value  of  the  current  density.  For 
a  given  bulk  concentration,  the  increase  in  the  applied  over- 
potential  increases  the  electrochemical  reaction  rate  on  the 
electrode  surface.  The  effect  of  the  applied  overpotential  on 
the  current  density  can  also  be  observed  in  Figures  3  through  5> 
In  each  of  these  figures,  a  plot  of  the  limiting  current  - 
density- vs .-time  curve  is  also  shown.  The.  limiting  current 
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zero  at  any  tine.  After  examining  Figures  3  through  5»  find 
that,  for  the  case  of  -0.30V,  some  of  the  values  of  current 
density  exceed  the  limiting  current  density.  This  phenomenon 
becomes  more  significant  when  the  bulk  concentration  of  pota¬ 
ssium  nitrate  is  decreased.  The  reason  for  this  phenomenon 
probably  is  that  at  the  high  applied  overpotential ,  many 
reactions  can  subsequently  be  expected  to  take  place.  The 
maximum  current  appearing  in  the  curve  for  the  case  of  -0.50V 
is  a  result  of  this  behavior. 

When  the  applied  overpotential  is  fixed,  the  hetero¬ 
geneous  reaction  rate  constants  in  both  the  cathodic  and  anodic 
directions,  K,  and  K2 ,  are  fixed.  The  supply  of  the  nitrate 
ions  on  the  electrode  surface  comes  only  in  the  way  of  diffusion 
of  the  nitrate  ions  from  the  bulk  solution.  Increasing  the 
bulk  concentration  will  increase  the  rate  of  nitrate  ion  supply 
to  the  electrode  surface,  thusly  increasing  the  current  density. 
Figures  7  through  9  show  the  effect  of  the  bulk  concentration 
of  the  electrolytic  solution  on  the  current  density  with  fixed 
applied  overpotential.  From  the  three  figures,  we  can  see  that 
the  current  density  is  approximately  proportional  to  the  bulk 
concentration  of  potassium  nitrate,  except  for  the  case  of  -C: . 

When  electrolytic  solution  contains  cadmium  ions,  a 
different  set  of  results  was  obtained.  The  experimental 
current-density-vs . -time  curve  obtained  in  a  J.C  02pl-l 
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Effect  of  the  Bulk  Concentration  of  the  Electrolytic  Solution 
on  the  Current  Density  with  Fixed  Applied  Overpotential . 
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Effect  of  the  Bulk  Concentration  of  the  Electrolytic  Solution 
on  the  Current  Densi  ty  with  Fixed  Applied  Overpo tential . 
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Cd(N0^)2  solution  and  at  -0.40V  applied  overpotential  is  shown 
in  Figure  6.  In  addition  to  the  electrochemical  reduction  of 
the  nitrate  ions  on  the  electrode  surface,  a  precipitation 
reaction  also  takes  place  in  the  solution  because  the  cadmium 
ions  in  the  electrolytic  solution  can  coprecipitate  with  the 
hydroxy  ions  produced  by  the  electrochemical  reaction  on  the 
electrode  surface.  The  current  density  is  expected  to  be 
greater  than  that  when  the  cadmium  ion  is  absent  from  the 
electrolytic  solution.  Figure  10  shows  the  current-density - 
vs. -time  curves  obtained  in  a  0.005WI  KNO^  solution  and  in  a 
0.0025M  Gd(NC^ ) 2  solution.  The  applied  overpotential  is  -0.2-07 
for  both  cases.  From  the  curves,  one  can  see  that  the  current 
density  obtained  in  the  0.0025M  Cd(N0^)2  solution  is  larger 
than  that  obtained  in  a  0.005M  KNC^  solution.  The  reason  for 
this  is  that  the  consumption  of  the  hydroxy  ions  due  to  the 
precipitation  reaction  establishes  a  driving  force  to  produce 
more  hydroxy  ions  and  thus  promote  the  electrochemical  reaction. 
Some  interesting  facts  may  be  seen  from  the  figure.  Firstly, 
the  difference  between  two  current  densities  increases  as  time 
increases,  finally  reaching  a  constant  difference;  secondly, 
the  rate  of  decrease  of  the  current  density  with  time  is  slower 
for  the  case  of  using  Cd(N0^)2  as  the  electrolytic  solution. 

One  may  thus  obtain  a  larger  steady  state  current  density  using 
C d ( N 0 ^ P  as  the  electrolytic  solution. 
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CHAPTER  5 

NUMERICAL  SOLUTIONS 
OF  THE  THEORETICAL  MODEL 

5 • 1  Orthogonal  Collocation  Method 

There  are  many  numerical  techniques  which  can  be  used 
to  solve  the  differential  equations  of  the  theoretical  model 
presented  earlier  in  Chapter  3* 

Orthogonal  collocation  method  was  selected  as  the  method 
of  solution.  The  orthogonal  collocation  method  requires  very 
little  implicit  mathematics  and  results  in  a  large  savings  in 
computer  time  over  the  common  finite  difference  scheme  such  as 
the  Crank-Nicholson  method.  A  brief  description  of  the  ortho¬ 
gonal  coxlocation  method  is  presented  in  Appendix  C. 

5*2  Working  Equations 


5*2.1  Dimensionless  Form  of  the  Differential  Equations 


The  orthogonal  collocation  method  requires  the  region  of 
solution  to  be  restricted  to  the  interval  C 3 , 1 3  .  This  is 
accomplished  by  normalizing  the  spatial  variable  by  the 
parameter ,  £ ,  the  diffusion  thicxr.ess,  ^  =x/ /  . 


dimensionless  forms* 

*  *  ■'-> 
r*  _  —  ,  n  —  - 

1  cl  2  c? 

1  1 

the  dimensionless  time  is  defined  as: 


Substituting  the  above  dimensionless  variables  into  equations 
(17)  and  (19)  and  rearranging  leads  to  the  following  dimension¬ 
less  eauations: 


.(.Eulh 
l  V^2 


*  2„  * 

^  ^  '‘l  ,  ,  _  *  „  *q  3?Ist 

•2T  ‘s'?2  '  r  1  T'“2  ;"777r 


where  K=  -i-£—  and  2Ksp=  -ASF  , 

32 

are  the  dimensionless  homogeneous  reaction  rate  constant  and 
solubility  product,  respectively.  Similarly,  equations  (21) 
and  (22)  become: 


U1  \  1>  U1 


*  ?  * 

z:2  _  yc2 


V>J 
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The  initial  conditions  are: 


1.  att*0,  C.  =1, 


for  all  0 


s  l- 


2.  atr« 


=0,  C2"=  C°/C°  ,  for  all  0  <  £  < 


1. 


The  first  two  boundary  conditions  come  from  the  dimensionless 
form  of  equations  (32)  and  ( 33 ) 1 


D1  C°  f  >Gi 
1.  (  1  0  1 


*  N 


i  J^o 


K2  (C°)3  <C°V  -  (KjC“)  cf' 


=  0  , 


31  f9c-  1 

2.  3(  — )  ^r?- 

C2  L23  J$>0 


f  3C; 


=  0 


:=o 


The  other  two  boundary  conditions  are: 


3 .  at  ^  =1 ,  C1  =1  , 


for  ^  >  0  , 

4.  aty«l,  C2*=C^/cr,  forf>0. 


5.2.2  Discretized  Equations 


The  solution  is  approximated  by  a  (L+2 ) th-order  Legendre 
polynomial  which  satisfies  the  differential  equations  and 
boundary  conditions  exactly  at  N+2  points.  Those  points  are 
chosen  to  be  zeros  of  the  Legendre  polynomial  of  degree  over 
the  interval  (0,1)  and  with  the  boundary  points  0  and  1. 
Discretizing  the  spatial  derivatives  at  those  points  leads  to 
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the  following  set  of  first-order  ordinary  differential  equations 


i=  2, . N  +  l.  (39) 


1—  2,#..#N+1,  (  irO  ) 

where  ^  ^ ,  i=2,...  N  +  l  are  the  N  inferior  collocation  points 
while  1  and  ^+2  are  exterior  collocaxion  points  which 
correspond  to  the  boundary  points  ^  =0  and  =1, 
respectively,  and  0=  ^  1  <  ^  ^  N+2  =  1.  i=2,  .  . .  N  +  l.  3j_  , 

is  a  (N+2)  by  (N+2)  coefficient  matrix  which  comes  from  the 
discretization  of  the  second  derivative  (  d2C1*/^'^2  or 
O  C2  / S  *5  )  at  each  collocation  point.  A  detailed  explanation 
of  matrix  3  is  shorn  in  Appendix  C.  C,  (^  .,*£■)  and  C~  (t  ,,t) 
are  the  concentrations  of  species  1  and  2,  respectively,  ac 
time  'C  and  position  ^  ^ . 
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where  ^  is  an  element  of  another  ( M -*-2 )  by  (N+2)  coefficient 
matrix  which  comes  from  the  discretization  of  the  first-order 
derivative  (  2>C^  or  "^C2  )  at  each  collocation  points . 

Since  only  boundary  conditions  1  and  2  (equations  (43)  and  (44)) 
require  the  evaluation  of  the  first-order  derivative  at  ^ 1 =0 , 
only  the  first  row  of  the  matrix  A,  i.e.,  A,  is  used  here. 

1  I  o 

The  initial  conditions  for  both  cases  are  the  same: 

I.C.  's 

1.  at  -f=0,  c1*(^i.T)=  1>  i  =  l,...N+2, 

2.  at  ^  =0,  c2*('|i,^)=  C^/C^  ,  i=l,..,N+2. 


5*2.3  Calculation  of  Current  Densit\ 


The  current  density  is  related  to  the  flux  of  the  nitrate 
ion  at  the  electrode  surface  by  equation  (34).  The  discretized 
and  dimensionless  form  of  equation  (34)  is; 
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3 . 3  The  Introduction  of  a  Solir.e  Point;  to  the  Tiscretited 
Equations 


The  diffusion  thickness  introduced  in  the  derivation 
of  the  dimensionless  form  in  the  previous  section  can  he 
justified  by  selecting  the  total  electrolysis  time  such 
that  diffusion  effects  are  negligible  at  x=£  during  chat 
period.  At  the  beginning  of  the  electrolysis,  the  diffusion 
thickness  selected  according  to  the  total  electrolysis  time  is 
too  large.  This  leads  to  a  concentration  profile  which  is 
fiat  in  most  of  the  interval  with  a  sharp  drop  in  a  very  small 
region  near  *5-  =0.  This  would  require  a  large  number  of 
collocation  points  leading  to  a  very  stiff  set  of  ordinary 
differential  equations. 

This  difficulty  can  be  overcome  by  the  use  of  a  method 

f  2  Q  \ 

called  spline  collocation'  '  based  on  the  diffusion  boundary 
concept.  This  method  maintains  a  fixed  low  number  of  collo¬ 
cation  points  in  the  interval  which  is  very  small  initially 
and  will  be  increased  as  time  increases  to  account  for  the 
thickening  of  the  diffusion  layer.  In  the  regions  outside 
this  interval,  the  concentration  is  flat.  As  a  matter  of 
fact,  the  spline  point  can  be  located  as  close  to  the  inter- 

(  ?  Q  ) 

face  as  to  achieve  any  desired  accuracy  , 

t  the  concentration  gradient 


In  order  to  be  sure  tha 


at  the  spline  point  is  zero,  the  concentration  at  the  Nth 
collocation  point,  i.e.,  the  last  one  before  the  spline 
point,  is  tested  to  assure  that  it  is  within  a  very  small 
range  of  the  bulk  condition.  At  the  time  when  such  a  test 
fails,  the  spline  point  is  moved  further  into  the  solution, 

Introducing  a  spline  point, ^ s ,  G<  ^ c k  1 ,  a  new  variable, 
z= 'V /  "j:s,  is  introduced.  The  discretized  equations  we  obtained 
for  the  case  of  solution  containing  cadmium  ions  are: 


dC 

1 

df  la, 


_  ^  r  :;+2  * 

*  <  v2  >  <  r*  >  21 


i=2>  «  e 


a  (  >2  ) 

<X+2  *  *-i 

T  3,-  ,C2  (a,.T) 

-2K  f 

z  •  5  s 

^  *  ■  ■  -  f  J  ^  J  J 

j  =  l 

l 

^C2  (zrX)) 


rV  o  T) 

*  2 


1=2  ,  . , .h*l . 


For  the  case  of  solution  containing  no  cadmium  ion, 


c  n 
-/  ( 


Equations  (43)  through  (46),  which  cone  from  the  four  c cur. car; 
conditions,  become: 

D  cb  ?l+2 

la  (~fs  }  (  )[2I  Aifjci  (sj»T  J]  +  [k2(cJ)3c2  (  =  1,^) 

-  (K-)c;(t1>T)]  =0,  for  allt>0.  '  (32) 


2.  3(  • 


Di  .  <  ?il2 


L>[Z  ai,3Ci#(z3^0  +  [E  ai,;c2*^*^ 

2  U  j  =  l  ^  c  -  =  l  J 


for  all  X>°* 


3'  ci*(^s.t)“  ci*(::;.2't)=  for  ail  T>o, 

4*  52*^a-T>*  s2#‘*a»a»t)*«I/«i.  aii  t?c 


The  initial  conditions  become: 
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A 


4 


T  r  '  <3  • 

1  •  ^  •  O  « 

* 

1.  at  ^=0,  C1  (z^T)3  1.  i=l,...N+2, 

2.  at  r=0,  C2*(ZjL,r)=  C2/C2,  1=1,  .  ..N+2. 


The  current  density  should  he  changed  to  s 


5.4  Solution  Procedures 


For  the  case  of  a  solution  containing  cadmium  ions, 
equations  (48),  (49)  and  (52)  through  ( 55 ) »  with  the  initial 
conditions,  are  the  working  equations  used  to  solve  for  the 
concentration  profiles  of  N'O^-  and  OH-  ions.  The  Cd"1"^  ion 
concentration  is  calculated  by  the  electroneutrality  conditic 
in  the  solution: 
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For  the  case  of  a  solution  containing  no  cadmium  ion, 
equations  (50),  (51)  and  (52)  through  (55)  with  the  initial 
conditions,  are  the  working  equations. 

Equation  (56)  is  used  to  calculate  the  current  density 
at  each  time  /f  for  both  cases. 

5 • 5  Computer  Program  Structure 


The  computer  program  for  solving  the  discretized 
equations  consists  of  a  main  program  and  5  main  subroutines. 
The  first  main  subroutine  evaluates  the  collocation  points 
(roots)  of  the  Jacobi  polynomial  of  order  N,  as  well  as  the 
firs  and  second  derivatives  jf  tne  polynomial  at  the  roots. 
These  derivatives  are  then  used  in  a  second  main  subroutine 
to  caxculate  the  discretization  coefficient  matrices  A  and  3. 
The  third  main  subroutine  is  used  to  perform  a  semi-implicit 


integration  using  Gear's  routine 


(29) 


Inis  subroutine  cans 


4  external  subroutines.  The  first  of  these  contains  the 

explicit  expressions  (48)  and  (49)  (or  (50)  and  (51)  in 

another  case)  for  the  discretized  coupled  first  order  differentia- 

equations  containing  the  3.  .  terms;  the  second  contains  the 

->  j 

Jacobian  matrix  for  the  same  equation;  the  third  performs  the 
first  stage  of  Gaussian  elimination^  ;  of  the  Jacobian  matrix; 
tr.e  _ast  performs  the  back  substitur.cn  of  Gaussian  elimination 


of 
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of  the  Jacobian  matrix.  The  first  external  subroutine  called 
by  the  integration  subroutine  also  calls  another  subroutine 
coming  from  the  IMSL  package.  This  subroutine  solves  equations 

(52)  and  (53)  simultaneously  to  give  the  surface  concentrations 

*  * 

of  C.  and  C-  .  The  IMSL  subroutine  calls  another  function 

1  c. 

which  feeds  the  equations  to  be  solved.  Before  writing  the 

*  * 

function,  C ^  (z^,^)  (surface  concentration  of  Cp  )  is  solved 

-g.  -fr 

in  terms  of  (z-,^),  i  =  l,.,.Ni-2,  and  (t-.'T),  i  =  l ,  .  .  .11+2  , 

using  equation  (53)  »  and  then  substitute  into  equation  (52) 

* 

to  eliminate  (z^,^).  Thus,  only  one  equation  which 

* 

contains  only  one  unknown,  is  needed.  Cnee  we 

*  * 

obtain  (z^, '-£-),  C?  (z^,^)  can  be  calculated  easily. 

The  Jacobian  matrices  for  both  cases  can  also  be  evalu- 

* 

ated  after  substituting  the  expression  of  (z.,,^)  into 
equations  (43)  and  (49)  (or  (50)  and  (51)  in  another  case). 

For  the  case  of  solution  containing  no  cadmium  ion,  the  Jacobian 
matrix  is: 
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which  is  a  2M  by  2N  matrix.  For  the  case  of  solution 
cadmium  ions,  the  Jacobian  matrix  is  the  same  except 
the  terms  in  the  diagonal  line  of  the  third  quadrant 
subtracted  by  X  and  the  terms  in  the  diagonal  line  of 
fourth  quadrant  are  all  subracted  by  2K( f*( DDXsp/C^  v 


containing 

that 


are  all 
the 


T) 


The  values  obtained  from  the  integration  subroutine  are 
then  used  in  a  fourth  main  subroutine  to  calculate  the  current 
density  by  equation  (56)  .  The  last  main  subroutine  is  used  to 
evaluate  the  concentration  at  any  desired  point  between  0  and  1 
by  interpolation  using  the  values  at  the  collocation  points. 

The  structure  of  the  whole  program  is  shown  in  Figure  11. 
Programs  from  both  cases  are  listed  in  Appendix  D,  Programs 
herein  were  executed  on  an  AMDAHL  4?0/V6  computer. 
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—j  i  tt;  X 
U  l  *r»i  4.— '-J 

F  ITT  I?*’ G  OF  THE  THEORETICAL 
MODEL  WITH  EXPERIMENTAL  DATA 

6 . 1  The  Determination  of  the  Electrochemical  Reaction 
Kinetics  Parameters 

The  unknown  heterogeneous  reaction  rate  cons¬ 
tants  wnicn  appear  in  equation  (52)  were  dexermined 
from  xhe  currenx-density-vs . -time  data  obtained  in  a 
solution  which  did  not  contain  cadmium  ion.  The  di¬ 
fferential  equations  (50)  and  (51)  which  describe  the 
transport  process  were  solved  by  computer  to  determine 
the  rate  constants.  Comparing  equation  (31)  with  equa 
tion  (31a),  one  can  see  that  both  K1  and  K9  depend  on 
the  applied  potential  with  exponential  dependence. 
Thus,  X.  and  ^  should  be  constant  when  the  applied 
potential  is  maintained  at  a  constant  value  during  the 
electrolysis.  K1  and  K^  were  obtained  by  comparing  the 
experimental  data  in  the  absence  of  cadmium  ion  in  the 
electrolyte  with  the  computer-calculated  current  der.si 

The  experimental  current-der.sicy-vs  , -time  data 


data 


with  the  electrolyte  containing  0.005'^  nitrate  ions 
and  the  applied  overpotential  of  -0,407  was  first 
fitted  to  the  computer-simulated  data  to  obtain  the 
K1  and  values.  Due  to  the  lack  of  information 
about  the  values  of  K1  and  ,  rough  estimates  of  X. 
and  K0  were  evaluated  by  substituting  the  applied 
potential  V,  which  is  equal  to  the  value  of  the  applied 
overpotential  plus  the  rest  potential,  to  equations 
(31)  and  (31a).  The  estimated  values  of  X,  and  K,  were 
increased  or  decreased  until  the  computer-calculated 
current-density-vs .-time  curve  had  the  same  shape  as 
the  same  curve  obtained  experimentally,  ('iote,  the 
background  water  decomposition  current  has  been  sub¬ 
tracted.  ) 

After  we  obtained  the  right  shape,  the  anodic 
reaction  rate  constant,  X^ ,  was  then  fixed.  By  changing 
X,  ,  one  can  make  the  entire  curve  shift  upward.  Con¬ 
versely,  decreasing  the  K,  value  will  decrease  the 
current  density  at  each  time  point  and  thus  shift  the 
entire  curve  downward.  A  set  of  current-density-vs . - 
time  curves  with  fixed  K7  values  and  different  values 
of  X.  are  plotted  in  Figure  12.  The  values  of  X,  and 
K0  when  the  concentration  of  nitrate  ions  is  0 . CC 
and  the  applied  overpotential  is  -G.^CV  w as  obtained 


Electrolysis  Time,  ms 
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by  interpolating  the  experimental  curve  between  those 
curves  in  Figure  12.  The  values  of  X1  and  X 2  which 
achieved  the  best  approximation  of  the  experimental 
data  are  K.  =  0.12  X  10~-^  and  K^=  0.6  xlO--5,  respectively 

These  values  were  then  used  to  predict  the  cu- 
rrent-density-vs . -time  curves  at  the  other  ion  con¬ 
centrations.  Figure  13  shows  the  computer-simulated 
current  density  curves  and  experimental  curves  at  va¬ 
rious  nitrate  ion  concentrations . 

The  values  of  and  K2  at  the  applied  overpo¬ 
tential  of  -0.60V  were  different  from  those  at  -0.40V. 
The  same  procedure  was  used  in  evaluating  the  values  of 
and  K2  at  the  condition  of  -0.60V.  Figure  14  shows  a 
set  of  calculated  current-density-vs . -time  curves  with 
fixed  K2  valued  and  different  values  of  K1 .  The  value  o 
K„  at  the  condition  of  -0,60V  was  then  obtained  in  the 

i. 

same  way  as  before.  The  values  of  K,  and  K0  were  deter- 
mined  to  be  0.10  x  10  J  and  0.50*10  ,  respectively. 

These  values  were  then  used  to  predict  the  current-den¬ 
sity-vs  .  -time  curves  at  the  other  ion  concentrations. 
Figure  15  shews  the  comparison  of  the  predicted  and  ex¬ 
perimental  current-density-vs , -time  curves. 

The  reaction  rate  constants  at  the  applied  overpo 
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Figure  I'Jt  Experimental  and  Computer-Simulated  Curren  t-Densi  ty-vs  .  -T  ime 

Curve:;  with  Various  Nitrate  Ion  Bulk  Concentrations  in  the  Case 
of  FI  t-c  I  roly te  Containing  No  Cadmium  Ion. 
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tential  of  -0.30V  were  not  pursued.  The  reason, which 
was  discussed  previously,  is  that  at  this  high  applied 
overpotential  more  complicate  sequences  of  reactions 
take  place.  This  simple  reaction  expression  can  not 
describe  the  phenomena. 

The  heterogeneous  rate  expression  so  determined 
is  believed  to  be  a  correct  one.  This  claim  is  supported 
by  the  fact  that  the  rate  constants,  and  Kg,  are 
independent  of  the  bulk  concentration  at  a  given  poten¬ 
tial.  This  is  indeed  the  case  as  was  shewn  in  Figures 
13  and  13. 

6.2  The  Determination  of  the  Homogeneous  Precipitation 
Reaction  Rate  Constant 

Higher  current  density  was  observed  when  cadmium 
nitrate  was  used  instead  of  potassium  nitrate  as  the 
electrolyte  in  the  electrochemical  cell.  The  hydroxy 
ion  produced  by  the  reduction  of  nitrate  ion  copreci¬ 
pitated  with  the  cadmium  ion  in  the  solution.  This 
increased  the  reduction  reaction  rate  and  led  to  higher 
current.  This  homogeneous  reaction  was  assumed  to  be 
linearly  dependent  or.  the  degree  of  supersaturation  of 
cadmium  hydroxide  (see  Equation  (1^)).  Since  the  beta- 


roger.eous  reaction  rate  constants  X„  and  were  deter¬ 
mined,  the  homogeneous  rate  constant  d  could  be  evalu¬ 
ated  'ey  fitting  the  second  set  of  experimental  data 
with  the  theoretical  model  to  determine  the  value  of 

Figure  16  shows  the  current-density-vs , -tine 
curves  with  various  d  values  while  the  other  parameters 
were  held  at  a  constant  value.  The  current-density-vs . - 
time  data  obtained  in  a  0.0025M  cadmium  nitrate  solution 
at  the  applied  overpotential  of  -0.407  is  also  plotted  in 
Figure  16.  The  reaction  rate  constant,  > ,  was  then  eva¬ 
luated  approximately  from  Figure  lc.  The  most  probable 
homogeneous  reaction  rate  constant  was  determined  to  be 
SO  sec-1. 


k=  100 


ity-vi; . -Time  Curves  with  Fixed 
ua  Reaction  Rate  Coikj tault;. 
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CHAPTER  7 

DISCUSSIONS  AND  CONCLUSIONS 

7 • 1  Surface  Concentrations  of  Various  Ions  as  Junctions 
of  Time 


The  knowledge  of  the  surface  ion  concentrations  can  pro 
vide  information  concerning  the  electrochemical  deposition 
process.  Although  the  surface  ion  concentrations  were  not 
available  from  the  experiment,  that  information  could  be 
generated  by  computer  simulations. 

Figure  17  shows  the  surface  concentrations  of  N0^~  and 
0H~  ions  as  functions  of  time  in  the  absence  of  cadmium  ion 
in  the  solution  at  two  potentials,  -O.lOV  and  -0.60V  from 
equilibrium  potential  and  at  a  N0^~  bulk  concentration  of 
0.005M.  In  generai,the  surface  concentration  of  NO,"  ions 
decreases  as  time  increases,  while  the  surface  concentration 
of  0H~  ions  increases  as  time  increases.  This  is  due  to  the 
electrochemical  reaction  (23),  in  which  one  nitrate  ion 
reacts  with  two  electrons  and  which  produces  three  hydroxy 
ions . 

For  the  case  of  -O.m-CV  overpotential,  the  CH-  concen¬ 
tration  increases  sharply  from  the  bulk  value  of  nearly  cere 


KEY 


?4 
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to  about  0.5  dimensionless  concentration  units  in  less  than  1 
msec  before  the  hydroxy  ions  generated  are  removed  effectively 
by  diffusion.  Thereafter,  the  OH-  concentration  increases  at 
a  slower  rate.  For  ions,  the  surface  concentration  drops 

to  about  0.72  dimensionless  concentration  units  from  the  bulk: 
value  of  1.0  in  less  than  1  msec  and  then  decreases  at  a 
slower  rate  when  the  diffusion  can  effectively  supply  the 
reactant  for  the  electrode  reaction.  The  sharp  changes  of  the 
surface  concentrations  of  and  CK  ions  cause  the  sharp 

decrease  of  the  current  density  in  a  short  time  period  in 
the  beginning  of  the  electrolysis. 


At  a  more  negative  overpotential  of  -C.60V,  the  increase 
in  surface  CH-  concentration  and  the  decrease  of  NO^  concen¬ 
tration  are  even  more  significant  and  last  over  a  longer  time 
period.  The  surface  concentration  of  NC^  drops  to  about 
0.5  dimensionless  concentration  units  ar.d  OH  concentration 
increases  to  about  G.c7  dimensionless  concentration  units  in 
2.0  msecs.  The  behavior  is  the  consequence  of  the  higher 
overall  reaction  rate  at  -0.60V.  The  nitrate  ions  were  con¬ 
sumed  at  a  higher  rate  which  produced  more  hydroxy  ions, 
longer  induction  time  is  needed  for  diffusion  to  effectively 
supply  the  reactant  from  the  bulk  solution  and  remove  tr.e 
oroduct  from  the  electrode  surface. 


Similar  profiles  of 


suriacs 


cor 


and  CH 
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at  lower  bulk  NO^~  concentrations  -were  also  obtained.  These 
are  shown  in  Figures  18  and  19* 

When  the  electrolyte  contains  cadmium  ions,  the  homo¬ 
geneous  precipitation  reaction  of  cadmium  ions  and  hydroxy 
ions  near  the  electrode  surface  promotes  the  electrochemical 
reaction  on  the  electrode  surface  in  the  cathodic  direction. 
This  causes  the  nitrate  ions  to  be  consumed  at  a  higher  rate 
than  the  case  when  no  cadmium  is  present.  This  is  shown  in 
Figure  20  where  the  homogeneous  reaction  rate  constant  is  50 

4 

sec~l.  It  is  interesting  to  note  that  the  cadmium  ion  surfac 
concentration  increases  initially  and  then  starts  to  decrease 
This  phenomenon  is  explained  as  follows.  Immediately  after 
the  electrolysis  began,  N0^~  ions  were  consumed  by  the  elec¬ 
trochemical  reaction  which  produced  an  N0^“  concentration 

gradient  made  the  N0^~  ions  in  the  bulk  solution  diffused 

+2 

toward  the  electrode  surface.  The  Gd  ions,  on  the  other 

hand,  moved  with  the  NO^-  ions  in  the  same  direction  in  order 

to  maintain  the  elecroneutrality  condition.  In  the  meantime, 

the  OH  concentration  was,  however,  not  so  high  as  to  consume 
+2 

all  the  Cd  ions  brought  in  by  the  transport  process.  The 
cadmium  ions  was  thus  accumulated  and  resulted  in  a  higher 
value  than  its  bulk  condition.  As  scon  as  the  surface  c oncer, 
tration  of  CH  ion  raised  to  some  value,  the  surface  concen¬ 
tration  of  the  Cd  ^  icr.s  starts  ~o  decrease  because  of  the 
abundant  supply  of  OH"  ions  consuming  the  cadmium  ions  at 
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this  time. 


For  the  nitrate  ions,  the  surface  concentration  also 
decreases  as  time  increases  at  a  rate  somewhat  larger  than 
that  for  the  case  of  solution  containing  no  cadmium  ion.  This 
is  shown  in  Figure  21.  Figure  21  also  shows  the  difference  in 
the  surface  concentration  of  the  hydroxy  ion  between  those  two 
cases,  i.e.,  in  the  presence  and  in  the  absence  of  cadmium 
ions.  The  initial  behaviors  of  two  cases  are  very  similar. 
After  this  induction  period,  the  ion  concentration  for  the 
case  of  a  solution  which  contains  cadmium  ions  decreases  at 
a  more  gradual  rate  as  a  result  of  the  homogeneous  precipi¬ 
tation  reaction  between  cadmium  and  hydroxy  ions. 

7-2  Concentration  Profiles  of  Various  Tons  at  Various  Time 


For  the  case  of  solution  containing  no  cadmium  ion,  the 
nitrate  ions  were  consumed  and  the  hydroxy  ions  were  produced 
on  the  electrode,  the  concentration  gradients  were  established 
which  made  the  nitrate  ions  diffused  from  the  bulk  solution 
toward  the  electrode  surface  and  the  hydroxy  ions  diffused 
away  from  the  electrode  surface  toward  the  bulk  solution. 
Figures  22  through  27  show  the  concentration  profiles  of 
and  OH  ions  with  different  applied  overpoter.tials  at  some 
selected  times.  From  these  figures,  one  can  see  that  the 
diffusion  thickness,  £  ,  increases  with  time.  The  diffusion 
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Figure  22  Concentration  Profiles  of  Nitrate  and  Hydroxy  Ions  at 
Various  Time  for  the  Case  of  Solution  Containing  No 
Cadmium  Ion. 
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Figure  24s  Concentration  Profiles  of  Nitrate  and  Hydroxy  Ions  a 
Various  Times  for  the  Case  of  Solution  Containing  No 
Cadmium  Ion. 
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Figure  26»  Concentration  Profiles  of  Nitrate  and  Hydroxy  Ions  at 
Various  Times  for  the  Case  of  Solution  Containing  No 
Cadmium  Ion. 
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centration  Profiles  of  Nitrate  and  Hydroxy  Ions  at 
ious  Times  for  the  Case  of  Solution  Containing 


thicknesses  at  various  time  are  listed  in  Tables  2  and  3> 

The  concentration  profiles  of  N0^“,  OH"  and  Ca ions 
at  various  time  for  the  case  of  solution  containing  cadmium 
ions  are  shown  in  Figures  28  and  29.  The  shape  of  NO  ~ 
concentration  profiles  is  the  same  as  that  for  the  case  of 
containing  no  cadmium  ion  but  the  N  0  ^  ~  concentration  on  the 
surface  is  somewhat  lower  than  that  for  the  case  of  contain¬ 
ing  no  cadmium  ion.  The  diffusion  thickness  increases  as  time 
increases , 

In  the  previous  discussion  of  the  surface  concentration, 
the  surface  concentration  of  0H~  ions  increases  in  the  beginn¬ 
ing.  During  that  time,  OK"  ions  diffuse  from  the  electrode 
surface  to  the  bulk  solution.  Figure  28  shows  -chat  the 
surface  concentration  increases  and  the  diffusion  layer  thick¬ 
ness  also  increases  as  time  increases.  Then  the  surface 
concentration  starts  to  drop  and  the  diffusion  layer  thickness 
stays  at  about  the  same  value.  This  is  shown  in  Figure  30. 
Table  4  snows  the  diffusion  layer  thickness  for  the  diffusion 
of  OH"  ions  as  a  function  of  time.  In  Table  4,  one  can 
see  that  the  diffusion  layer  thickness  stops  increasing  as  a 
result  of  the  decreasing  surface  concentration  of  Oh’"  ion. 

+ 2 

The  Cd  concentration  profiles  are  shown  ir.  Figure  29  > 
As  the  profiles  indicate,  for  each  specific  time,  when  the 
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Table  2:  Diffusion  Thickness  as  a  Function  of  Time 
for  the  Case  of  Solution  Containing  No 
Cadmium  Ion.  The  Applied  Cverpoten'ial  is 
-0.40V.  The  Nitrate-Ion  Bulk  Concentration 
is  Q.005M. 


Time  (sec) 


0.1235  X  10 


-5 


0.6121  X 10 


-5 


0.1758  X10 


-4 


0.3572  *10 


-4 


0.1033  XlO~3 


0.4737  X10"3 


0.1258  X10 


-2 


0.7417  XlO' 


0.3946  XlO 


-1 


0.7032  *10 


-1 


Diffusion  Thickness( cm) 


0.25  XlO 


-4 


0.45  X10"4 


0.80  X 10 


— 


0.15  x 10 


0.25  X 10 


-3 


0.50  x 10"3 


0.80  XlO' 


0.20  X10 


1  n-2 


0.45  X 10 


-2 


0.60  X 10 


-2 
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Table  3:  Diffusion  Thickness  as  a  Function  of  Time 
for  the  Case  of  Solution  Containing  No 
Cadmium  Ion.  The  Applied  Overpotential  Is 
-0.60V.  The  Nitrate-Ion  Bulk  Concentration 
Is  0.005M. 


Time  (sec) 


0.1313  xio"-5 


0.4982  X 10 


-J5 


0.1642  XIO 


rzr 


0.3318  xio" 


0.1298  xio"3 


0.4407  XIO"3 


0.1258  XIO 


-2 


0.8718  XIO 


-2 


0.3229  XIO 


-1 


0.7004  XIO" 


Diffusion  Thickness  (cm) 


0.25X10 


-4 


0,40  X 10 


0.70  X  lo-i; 


0.10  X 10" 


0.25  xio" 


0.45  XIO"3 


0.75  xio"3 


0.25  Xio 


-2 


0.45  XIO' 


0.65  Xio 


-2 
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30 i  Concentration  Profiles  of  Hydroxy  Ion  When  the  Surface 
Concentration  of  Hydroxy  Ion  Starts  to  Drop. 
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Table  4;  Diffusion  Thickness  of  0K~  Ion  as  a  Function 
of  Time  for  the  Case  of  Solution  Containing 
Cadmium  Ions.  The  Applied  Overpotential  Is 
-0.40V.  The  Nitrate-Ion  Bulk  Concentration  Is 
0.005M.  The  Cadmium-Ion  Bulk  Concentration  Is 
0.0025M. 


Time  (sec) 

Diffusion  Thickness ( cm) 

0.1236  XIO'3 

0.25  X10~4 

0.6128  X10"5 

0.45  X10~4 

0.1798  XIO-4 

0.80  XIO-4 

0.3513  xio"4 

0.15  xio'3 

0.1017  X io~3 

0.25  XIO-3 

0.4856  X 10'3 

0.50  X 10'3 

0.1148  XIO-2 

0.70  XIO-3 

1  0.8811  XIO-2 

0.85  X10~J 

0.4 050  X 10"1 

0.85  xio'3 

0.7031  XlO-1 

0.35  X 10'3 
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concentration  near  the  electrode  is  greater  than  the  bulk 

+2 

concentration  of  Cd  ,  its  value  versus  distance  from  the 
electrode  decreases  slowly  in  the  beginning,  then  faster,  and 
more  slowly  again  before  finally  reaching  its  bulk  value.  At 
the  time  when  the  concentration  values  near  the  electrode 
surface  are  still  larger  than  the  bulk  concentration,  the 
concentration  decreases  to  some  minimum  value  at  some  distance 
from  the  electrode  surface  and  then  increases  until  it  reaches 
the  bulk  value.  The  location  of  the  minimum  concentration 
point  moves  toward  the  electrode  surface  as  time  increases. 

At  sufficiently  long  time,  the  minimum  point  finally  reaches 
the  electrode  surface. 

7 . 3  Conclusions 


The  transport  model  for  the  electrode  kinetics  and  the 
homogeneous  reaction  presented  previously  has  already  explained 
very  successfully  the  characterestics  of  the  current  density  as 
a  function  of  time  for  both  the  case  of  solution  containing 
cadmium  ion  and  that  without  cadmium  ion.  This  model  was  then 
used  to  simulate  the  concentration  profiles  for  various  ions 
involved  in  the  reactions.  The  behavior  of  the  various  ions  on 
the  electrode  surface  and  in  the  solution  was  obtained.  The  ii 
fusion  processes  were  also  known  after  examining  these  concen¬ 
tration  profiles  carefully.  These  can  be  summarised  as  follows 


1.  The  current  density  decays  with  time  after  a  double-layer 
charging.  The  time  for  this  charging  is  so  short  that  we 
can  not  observe  it  cy  the  experiment.  The  decay  for  she 
current  density  is  a  result  of  the  depletion  of  the  nitrate 
ions  on  the  electrode  surface. 

2.  The  current  density  for  the  case  of  a  solution  containing 
cadmium  ions  is  greater  than  that  for  the  case  of  solution 
containing  no  cadmium  ion.  This  is  because  the  precipitaticr 
of  Cd(CH) 2  promotes  the  electrochemical  reaction  or.  the  elect 
rode  surface  by  consuming  OH”  generated  in  this  reaction. 

3.  The  current  density  depends  on  the  bulk  concentration  of 


nitrate  ions  and  the  applied  overpotential.  Increasing 
the  NO^-  bulk  concentration  or  the  applied  overpotential 
will  increase  the  current  density. 

k.  The  applied  overpotential  is  one  of  the  factors  which 

changes  the  reaction  rate  of  the  electrochemical  reaction 
on  the  electrode  surface.  Increasing  the  applied  over¬ 
potential  will  increase  the  overall  reaction  rate. 

5.  For  the  case  of  solution  containing  no  cadmium  ion,  the 
surface  concentration  of  N 0 ^ "  decreases  as  time  increases 
while  the  surface  concentration  of  OH"  increases  as  time 
increases.  This  is  because  NC^“  ions  are  consumed  to  produce 
the  OH-  ions  during  the  electrochemical  reaction  on  the 
electrode  surface. 

6.  For  the  case  of  solution  containing  cadmium  ion,  the 

surface  concentration  of  icn  also  decreases.  The 

decreasing  rate  is,  however,  faster  than  that  for  the  c--- 


9? 


of  solution  containing  no  cadmium  ion.  The  surface 
concentration  of  OH  ion  increases  in  the  beginning  arc 
reaches  a  maximum.  When  the  production  rate  of  OH-  due  to 
the  electrochemical  reaction  is  less  than  the  consumption 
rate  of  0H“  due  to  the  precipitation  reaction,  the  surface 
concentration  of  CH_  starts  to  decrease.  To  maintain  elec- 
troneutrality ,  the  surface  concentration  of  Cdr^  ions  will 
exceed  its  bulk  value  and  then  decrease  to  the  value  lower 
than  its  bulk  value. 

7.  The  diffusion  processes  for  the  case  of  solution  containing 
no  cadmium  ion  is  that  the  nitrate  ions  diffuse  from  the 
bulk  solution  toward  the  electrode  surface  while  the 
hydroxy  ions  diffuse  in  the  opposite  direction. 

8.  For  the  case  of  solution  containing  cadmium  ions,  the 
diffusion  directions  for  both  the  nitrate  ions  and  hydroxy 
ions  are  the  same  as  that  for  the  case  of  solution  containing 
no  cadmium  ion.  However,  the  cadmium  ions  diffuse  toward 
the  bulk  solution  in  the  beginning  and  then  finally  change 
direction  toward  the  electrode  surface. 

7  »4  Recommendation  far  Future  Works 

Although  the  model  has  successfully  explained  the  cases 
of  -0.407  and  -G.cQV,  iz  is  not  likely  valid  for  the  case  cf 
-G.oCY  due  to  the  unexpected  characterestics  of  the  ourrer.t- 
de.nsity-vs . -time  curves  and  the  comparison  of  the  current 


density  with  the  limiting  case.  Vftien  the  applied  cverpoter.tial 
is  as  high  as  -0.80V,  the  reaction  mechanisms  and  electrode 
behavior  should  be  identified  by  doing  some  further  experiments 
Furthermore,  the  current-density-vs . -time  data  at  -C.60V  applie 
overpotential  for  the  case  of  solution  containing  cadmium  ions 
is  not  available  yet.  It  can  be  obtained  by  repeating  the 
experimental  work  we  decribed  before.  Once  the  data  is  obtained 
the  model  can  be  used  for  comparison  with  the  experimental 
data. 


Finally,  the  configuration  of  the  electrode  in  the 
impregnation  process  is  not  as  simple  as  we  have  used  in  derivi 
the  model.  It  is,  in  fact,  a  porous  electrode  in  the  real 
case.  If  a  micro-porous  electrode  is  possible  to  obtain, 
the  experiment  can  be  repeated  to  obtain  the  current-density- 
vs. -time  data.  Cnee  the  data  are  available,  the  transport 
processes  can  be  examined  for  the  condition  of  porous  electrode 
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APPENDIX  A.  1 


THE  "BASIC"  MAIN  PROGRAM  USED  IN  EXPERIMENTAL  WORKS 
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10 

R  EM* 

k*s  Me.*:*:*:*:*:*:  Kite:* 

.  4<. :  tN  f  •  ■  fM  >c :  r: '.f : .  f::f; :  t'.  . 

1 3 

PEN* 

13 

REM* 

20 

REN* 

THIS  IS  THE  BA 

S I C  PROGRAM 

—  J 

REM* 

T I ME  TRANS I ENT  E 

MF'ER:  IMENT . 

30 

REM* 

26-0  IS  A  SAMPLING  SUBPOUT I N 

33 

REM* 

LANGUAGE •  WHEN  " 

MUG"  IS  EXE 

40 

PEM* 

SAMPLE  THE  CURRENT  OUTPUT  A 

-  cr 
"tj 

REM* 

INTERNAL.  THE  ANALOG  SIGNAL 

30 

PEM* 

INTO  A  DIGITAL  I 

ATA  BY  A  12 

33 

REM* 

IN  THE  MEMORY.  V 

HEN  THE  NI.-N 

60 

REM* 

PEACHES  A  PRESET 

UALUE THE 

64 

REM* 

TERMINATED.  THEE 

E  DATA  ARE 

6  6 

PEM* 

DECIMAL  NUMBERS 

Af  iD  FR I  NTED 

i-.X 

REM* 

THE  CAR  I MEL 2 S 

used  in  th i 

7"0 

REM* 

71 

REM* 

LOCK  S 

ETT I NG 

c  4 

REM* 

ID  DESAMPLING  TIME  I NT- 

“*^T 

REM* 

L,'=i.nnr-iriEL 

NUMBER  Ur  IE 

r 

REM* 

SAMP  LEZ- 

FROM 

60 

REM  * 

M  =NUME ER 

07  POINTS  T 

•"2 

REM* 

THE  ARGUMENTS 

IN  "MUG"  SU 

54 

REM* 

•56 

REM* 

I F  S=0  THEN 

EXECUTE  TEE 

•Z’Zf 

REM* 

IF  S=1  THEN 

GO  TO  MUGGE 

'50 

REM* 

N=NUMBER  OF 

UALU2E  PA EE 

92 

REM* 

M=NUMBER  OF 

UALUES  TO  £ 

H  -• 

REM* 

XOO'J  IE  THE 

INPUT  ARRAY 

4* 

<• 

f; 

M 

t. 

f: 
4< 
f: 
4- 
■  f 
4> 
+ 
4 
f 
+ 
t 
.+ 


PROGEfiM 


'■ 0 .-1  15  THE  INPUT  nrPiiV 

V's*J.>  15  THE  OUTPUT  CU7PEN7  FiRPriV 


*54  PEj*1;+; 

58  PEM:* 

57  REM* 

•58  PEM* 

*H'?  i^EM******^*****^’**  k.r.  ......  k.r  v,  ri- 

1 70  COfr  4r  I G 1  5CRES  4  ' 

18*3  DIM  'A '•  10!)  .•  V1’!  w00*0 !.* ..  Z ■•! 3000  ••  .•  C '•!  500  ■  .■  T 1  300  ' 

1 ?0  5-0 

200  N=3 

201  rEM* 


202 

PEM* 

204 

PEM*****  INPUT 

Z-  nTn  r..{Nv..r-.v 

203 

PEM* 

206 

PEM* 

2 1  -0 

PRINT  " CLOCK  EET . 

‘  .■  "enrlPLINL  I.IL  .N’T  E.- ;  :.Te" 

200! 

INPUT  NX.  O'  >  •  1  \  1  > 

240 

-  -  IN  "-.HANNEL  NUT' 

■iBER  •  "  •  "NUM2ER  07  f  0-I:<TS-  "  'E 

230 

INPUT  I  D  2.>  ■  M..  F 

231 

PEM* 

02 


.“'EM* 

234 

REM***-+ 

CFiLL  7HE  MUG 

SUBROUTINE  i^r-iD  STnT  i  TO  2ni‘;F’L2  «■»* 

233 

REM* 

230 

REM-*1 

200 

00 1 7 '.ML 

G S  >  N  .•  X  '■  0  :>  .■  M  .•  V 

;  ©  >  ;> 

201 

REM* 

i 

2.*Z‘2l 

REM* 

1 

! 

26^ 

REM***:* 

CGNUERT  7HE  Cl 

.■RREH7  Dfi7ri  I NT 0  DECIMhL  NUMBERS  (•••+■•*■•*••  i 

263 

PEM* 

2*o 

REM* 

1 

270  FOP-'  1=07 CM—  1 

200  I  I  >  = ...  V  <  I  >  -  i  5  >  -  1 6-'2043 

200  next  i 


300 

FEN  ****  C0HMEC7  7HE 

CGMPU7EP 

70  THE  LINE 

I N7EF: 

302 

REM****  DIVIDED  THE 

CURRENT 

3 ' i-'  7rl  E  L-UP  F  rt  C  2 

r  >r  c_n 

v  t-  t-l 

304 

"EM***  *  0.  0*3  CM2  7C 

EECCFjE 

i'HE  CURRENT  Z  EN 

Z  *  l  »’ 

i 

303 

PEM****  PR  IH7  CURREr 

7  DENS I 7 

V  US  7 1  ME  DFiTri 

i-i. 

306 

REM* 

307 

REM 

310 

CONFIG  C  77'.'  > 

320 

PR I NT  TnB <  1 3 >  ..  "71  ME "  • 

7  7i£i  '  .  43  j  .• 

"CURRENT  DEi'iS 1 7 

t  -  it 

r 

330 

PRINT  7 Fit-  C  1 6 >  ''  <  SEC  >  ' 

..  7  no  <  30  > 

.•  "  <  FiMP  CM2  >  " 

340 

PRINT JFRIN7 

330 

FOR  1=1  70  i ? 

3*0 

7' r.'=u  1000000*1 

370 

C1  I ,« =-2‘: .  I“1  .*  *~  0.1-  0. 

*3T3 

300 

rF  IN  1  1  FlB*.  iZ>  1  ..  1  •  i  ...  .  1 

s  50  :> ..  c 

I  > 

3  ?® 

NEXT  I 

400 

FOP  1=1  70  M-  20-  1 

410 

7'  I 1 0-C'OC 00  <  CO*  I ' 

420 

U1  1  = —  1  —UJ  *•  4.  i  '  <M_  LJ  .  2 

■jTTTT 

430 

FF  I  N7  7 F2  13  .  7 I  •  77 

l*  •  TU .  -■  C 

X 

440 

NEC:  7  I 

430 

END 

i 
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APPENDIX  A. 2 

THE  ASSEMBLY  SUBROUTINE  "MUG"  USED  IN  THE  MAIN  PROGRAM 


Adilvuss  Contents  Instructions  Comments 
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ov-v 


Contents  Instructions  Coiiiinent: 


0 j lx)  C9  RCT  Return  to  main  program 


APPENDIX  B 

LISTING  OF  EXPERIMENTAL  DATA 
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:asie 


3.3  Curran* -Density- vs .  -Tine  Tata  for  the  Case  of 
Solution  Containing  '<ro  Cadmium  Ion.  .titrate- 
Ion  3ulk  Concentration*  O.QCfM,  Initial  pH=  3.0, 
Aoplied  Overpotential*  -0.30V,  Rest  Potential 
=* -0.35^V. 


OC.S 
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Table  3. 


''•urrer.t-ij  er.  s  i  ty  -  vs 


Solution  Containing  ': 


-me  -a:a  lor  tr.e  -a 
o  Cadmium  Ion.  :iitr 
Ion  Bulk  Concentration5  0.C025M*  Initial  pH= 3.0, 
Applied  0verootential=-0 . 6C7,  Rest  Potential 

=*-0.35^* 


CiJ 
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m  I  ’  Pi  P 
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■a-ie  3.9  Current-2er.sity-ys  ,-Tiae  Bata  far  tr.e  rase  of 
Solution  Containing  N'o  Cadniun  Ion.  j’itrat 
Ion  Bulk  Concentration3  C .CC125M,  Initial 
Applied  Overuote.ntiai3  -O.cCV,  Best  Potential. 
=  -0.354V. 


CUPPEHT  I'E’r  15 1 T' 
•  nf-’P  CMC 


.  000 1 
.  0002 
.  0003 
.  0004 
.  0003 
.  0006 


.  OOOi  7 


.  OiOOiO 
.  0009 


104.043 
141.976 
1—3.  r?6 
1 1 1 . 1 77 
102. 333 


.  004 

33 . 3 1 1 2 

.  006 

3oi .  4433 

.  003 

23.3164 

.01 

22. 3116 

.012 

4 1 . 0336 

.014 

13. 1336 

.016 

13.4043 

.  0 1 3 

17.2773 

.  02 

16.3263 

•  0  7,  •=;  T- 

I 

i 


•*ii  fli 
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Table  3.12  Background  Current  Density  of  Water.  Initial 
pK=  3.0,  Applied  Overpotential3  -0,90V,  Rest 
Potential3  -0.35^. 


.  0007 

47. 7013 

.  OOOS 

43.9433 

.  0009 

41.316 

.  00 1 

33.6363 

.0011 

37 . 1 344 

.0012 

.  00 1 3 

34 . 1 796 

.  00 1 4 

33.4234 

.0013 

32.3016 

.  OO 1 3 

31.3504 

.  00 1 7 

30 . 7992 

.  00 1 3 

30 . 043 

.  OO  1 9 

29. 6724 

.  002 

23.9212 

.0v34 

22.536 

.  006 

19. 9063 

.  003 

13.0233 

.01* 

16.5263 

.012 

15.3996 

.014 

14.6433 

.016 

13.3972 

.013 

13.5216 

.  02 

12.7703 

.  022 

12. 3947 

.  024 

11.  c.436 

.  026 

1 1 . 6436 

.  023 

10.3923 

.  03 

1 0 . 3923 

.  Uw2 

10.5167 

.  034 

10.5167 

.  036 

9. 76561 

.  033 

9. 76561 

.04 

9.  39 

cl  O 
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Table  5.1c  Currert-Dersity-v; 


-  _  me 


z  o: 


L3e  or 


Solution  Containing  Cadmium  Ions.  Nirrare- 
lon  Bulk  Concencration=  0 .C023Mi _ Cadni 
Bulk  Concentration5  0.00125M*  Initial 

Applied  Overpotential  =  -Q.4CV,  Has*  Po - 

=  -0.35^V. 


■JS  P,  P 


O 
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Table  3.19  Current -Density- vs . -T 

Solution  Containing  Cadmium 


ata  tor  mhe  Case  ci 
ns.  Nitrate- 
Ion  Bulk  Concentration*  0 . 00125^1  Cadmium-Icr. 
Bulk  Concentration*  0.C00o25’fl»  Initial  pH=  3*Q» 
Apolied  Overcotential*  -O.^OV,  Rest  potential 
=‘-0.3  5^  V. 


t-i 


Bulk  C  on  c  er.  tr  a  t  i  on =  0.006251*1,  Initial  prl= 
Aopiied  Overooter.tial=  -0.807,  Rest  Potential 
=* -0.354V. 


M 
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The  method  of  weighted  residuals,  KWR,  is  a  general 
method  for  obtaining  solutions  to  equations  of  change,  in 
our  case,  Fick's  second  law.  In  the  MWR ,  one  assumes  a  trial 
function,  usually  a  set  of  weighted  polynomials , substitutes 
this  trial  function  into  the  differential  equation  and  then 
selects  the  coefficients  of  the  polynomial  terms  by  specify¬ 
ing  that  the  residual  be  zero,  on  the  average,  at  certain 
points.  If  one  evaluates  the  differential  equation  at  the 
zeros  of  an  orthogonal  polynomial ,  the  residual  will  of 
necessity  be  exactly  zero  at  these  collocation  points.  3y 
increasing  the  number  of  collocation  points,  the  trial 
function  would  satisfy  the  differential  equation  at  more  and 
more  points^ 1 ^ . 


(1)  3. A.  Finlayson  £  Richard  Bellman  (Id,),  'mathematics 

in  Science  and  Engineering',  Vol.  8?,  Academic  Press, 
New  York,  1972,  pp  9 . 
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with  the  boundary  conditions  of  the  form 

y(o,t)=g(t)  (3) 

y ( 1 , t  )=h( t )  (4) 

and  initial  condition  such  as 

y(x,0 )=p(x)  (5) 

For  a  mass  transfer  problem  D  would  be  a  dimension¬ 
less  diffusion  coefficient  and  f(y)  would  be  homogeneous 
reaction  term.  Geometric  considerations  indicate  that 

o£  =  £3=0  and  mathematical  considerations  indicate  that  a 

( 2 ) 

suitable  but  not  unique  trial  function  isv  ' 

N 

y(x,t)»(l-x)y(0,t)+xy(l,t)+x(l-x)ZIai(t)P.  .(x)  (6) 

1 

The  ? 4  1  of  equation  (6)  refers  to  one  member  of  a 
0  **  x 

complete  set  of  polynomials.  Table  C.l  lists  some  Legendr 
polynomials ,  the  special  case  of  equation  (1)  -where  ct  =  3 
and  their  roots  .  Note  that  a  Legendre  polynomial  of  degr 
has  N  real,  distinct  roots  or.  the  interval  0<x<  1. 


z.k.  Finlay son  -  Richard  Bellman  (Ed 
in  Science  and  Engineering' ,  Vol,  37, 
New  York,  1971,  pp  105 ■ 


'Mathematics 
idenic  Press, 
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TA3LE  C.  1 


Legendre  Polynomial  and  Their 


KOOtS 


(3) 


fl 

PN  - 1 

Roots  of  P.. 

iN  ■ 

1 

-l+2x 

0.500000000 

2 

1 -6x+6x2 

0.211324865 

0.788675135 

3 

-l+12x-30x2+20x3 

0.112701665 

0.500000000 

0.887298335 

4 

1-20x+90x2-140x^+70x^ 

0.069431844 

0.330009478 

0.669990522 

0.930568156 

(3)  John  Villadsen,  'Selected  Approximation  Method  for 
Chemical  Engineering  Problems',  Reproset,  Copenhagen, 
Denmark,  1970,  pp  A3  and  A4. 
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I x  is  evident  that  regardless  of  the  value  of  N  cho¬ 
sen  one  can  always  rewrite  eqn.(6)  in  the  general  form 

M  +  2 

y(x, t)»y  xJ"1b,(t)  (?) 

*-  '  ~  O 

<j  =  l 

One  could  at  this  point,  as  is  normally  the  case  in 

using  MWR,  obtain  the  coefficients  b-(t)  by  substitution  of 

i 

eqn.(?)  in  the  differential  equation  whose  solution  is  de¬ 
sired.  However,  the  computer  programming  is  greatly  simpli¬ 
fied  when  they  are  written  in  terms  of  the  solution  of  the 
differential  equation  at  the  collocation  points,  y(x. ,t) 
where  i=i , 2 , 3 • • .N+2 ,  and  each  x,  is  one  of  the  N  roots  of 
the  polynomial  or  one  of  the  two  boundary  values.  So  in 
terms  of  the  N+2  collocation  points  we  may  write: 

y(x, ,t)=  T’xP'1c  .(t)  (8) 

u.  ,  u 

iie  now  have  an  approximating  polynomial  function  which 
we  will  force  to  fit  our  partial  differential  eqn .(2)  at 
the  N+2  collocation  points  (i.e.  cnoose  the  b,(t)  so  that 
our  partial  differential  equation  is  satisfied).  In  order 
to  do  this  we  will  need  both  the  first  and  second  derivatives 
witn  respect  to  x.  These  derivatives  can  be  calculated  ex¬ 
plicitly  in  terms  of  x  since  we  have  already  defined  the 
dependence  of  y  upon  the  spatial  coordinate , eqn .( o ; and  (8). 

It  is  easily  seen  that  if  one  considers  ail  the  co- 


1 3^ 


liocation  points,  x^,  equation  (5)  represents  one  center  o 
a  set  of  equations: 


N+2 

y(x1 ,  i)=  y~ xj~1'o  j(  t)  (9) 

]=' 


Cne  can  see  that  the  left-hand  side  is  a  vector  of 

the  trial  solution  evaluated  at  the  collocation  points. 

•  « 

Since  the  term  involving  xv'-*  depends  on  doth  i  and  a 
matrix  will  represent  this  term.  The  b.(t)  will  naturally 

O 

be  a  vector.  If  vie  let  a  single  car  represent  a  vector  and 
a  double  bar  a  matrix,  then  the  set  of  earns. (9)  through;!! 
can  be  'written  as  one  equation  in  this  simplified  notation 

y ( t / =Q • b( t )  lily 

wnere  Q..  ;=xi‘~^ 

•  u 

The  first  and  second  derivatives  with  respect  to  x 
may  then  be  obtained  from  equation  ill): 
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— 

dx“ 


(1^) 


where  C 


dx “ ” 1  i  ,  n  d2xJ_1  I 

,d  d \*i  °  jx. 


it  is  now  si.T.pie  to  soive  lor  Ditj  tne  vector  ox  o.(t; 

j 

coefficients,  since  simple  rules  of  matrix  algebra  apply  in 
this  notation.  Equation  (12)  yields 


o  ( t )  =  Q  "  y  ( t ) 


Substituting  this  into  eqr.s.(13)  and (14)  gives 


dy(_ 


dx 


Q'x  y(t)=  A  y(t) 


( 15) 


*1  y  (  v )  _  -\  ,-i  ~  1  ,,  (  j.  \  _  p  rrrzr. 
dxc 


(16) 


where  A=  C  and  3=  D  Q-^  or  in  terms  of  the  i-th  collo¬ 
cation  point,  after  some  simple  matrix  multiplication, 


dy ( x , t )  _ 


MtZ 


dx 


=  y  :  i  y  V  x  •  ,  t ) 


i  J 


,2.  - . 

—  2 ^  1  j  J,  =  T~ 3  --,-c 

.4  >  “i,J 


dx 


J=l 
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One  can  therefore  express  the  spatial  derivatives  in 
terms  of  the  trial  function  at  the  N+2  collocation  points  . 

Equation  (16)  can  be  directly  applied  tc  some  diffu¬ 
sion  process  such  as  that  described  by  ecns.(2)  through  (5). 
Equation  (2)  becomes 


dy  ( t )  _ 
dt 


f(y( t) ) 


(19) 


where  the  vector  dy(t)/dt  represents  the  set  of  time  deri¬ 
vatives  at  the  N't 2  collocation  points,  'tie  note  that  there 
is  no  explicit  spatial  (x)  dependence  in  eqn.(l?).  3y  assum¬ 
ing  an  explicit  polynomial  in  x,  eqn.(6),  and  evaluating  it 
at  specific  values,  the  collocation  points,  the  problem  is 
reduced  from  a  partial  differential  equation  with  two  inde¬ 
pendent  variables  to  a  set  of  ordinary  differential  equa¬ 
tions  , 


One  can  expand  eqn.(19)  into  its  separate  members  to 
solve  the  problem.  For  the  i-th  collocation  point  we  can 
write,  according  to  eqr..(19), 


for 


l.C.3.  •  • 
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One  now  considers  the  conditions  at  the  boundaries 
x^=0  and  x^^l.  Since  these  values  are  known  from  the 
original  problem  i.e.  eqns.(3)  and  (4)  one  can  subtituxe 
the  boundary  conditions  in  eqn.(20)  to  obtain 


d7dt’  ~  =  *l3i.r*(t,+  3i,N-2*h(t)' 


I3i,jy(xj»t)]+  f(y(xi’t)} 

Ja| 


for  i  =  2,3,4. ,.N+1. 


One  should  note  that  g(t)  and  h(t)  are  included  so 
that  the  boundary  conditions  may  be  time  dependent  in  which 
case  g(t)  and  h(t)  would  be  known  or  easily  obtained.  The 
diffusion  problem  has  been  transformed  into  a  sex  of  N 
simultaneous  first  order  ordinary  differential  equations 
in  N  unknowns  with  initial  conditions  and  can  be  solved 
numerically  by  a  number  of  standard  methods  on  a  compuxer. 

Evaluation  of  the  current  at  an  electrode,  i=  nFAD|^-!  , 

Sx  jx=0 

can  be  obtained  directly  from  eqn.(17)  once  the  concentra¬ 
tion  profile  is  obtained. 
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PROGRAMS  USING  ORTHOGONAL  COLLOCATION 


METKC I 


1 
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J  **********************  ***************«*«**«******4****«*-**««*** 

¥  'J 

C 

c 

THIS  PROGRAM  SOLVES  EQUATIONS  (50)  THROUGH  (55)  n'li.i  I  »o 

X 

J 

INITIAL  COnDI  YIGN3  USING  SPLINE  COLLOCATION  M  ET  h  0  D  FG  A  T  :i  i 

c 

c 

CASE  OF  THE  ELECTROL YT  E  CONTAINING  NO  CADMIUM  ION. 

'Z 

-  c 

EQUATIONS  (DO)  AND  (51)  A  HE  DISCRZTI LED  FIRST -ORDER 

C 

c 

ORDINARY  DIFFERENTIAL  EQUATIONS  «HICH  ARE  SOLVED  BY  AN 

w 

c 

SEMI-IMPLICIT  INTEGRATION  SUBROUTINE  CALLED  DIFSUB  USING 

c 

c 

GEAR'S  ROUTINE.  EQUATIONS  (52)  THROUGH  (55)  ARE  SOLVED 

'w 

SI  J  ULI  AN  ED  US  L  Y  3  Y  A  SUBROUTINE  CALLED  l  SY3TM  Nil  1CH  COMES 

c 

V- 

FROM  I  MS  4>  SUBROUTINE  PACKAGE  u  OBTAIN  THE  SURFACE 

c 

c 

CD  NC  EN  IRA 1 .0 no  OF  NITRATE  AND  HYDRO*!  ION'S.  EQUATIONS  (54) 

c 

c 

AND  (5  5)  ASi  In  SIR  B  UL  K  VALUES.  TEE  CO.NCENTIi  ATION  PROFILES 

J 

c 

THUG  03IA1NLD  ARE  THEN  USED  TO  CALCULATE  THE  CURRENT 

c 

c 

DENSITY  jY  :t  UA  TIG  N  (5b).  TiiE  CURRE..T  DENSITY  DATA  ARE 

G. 

c. 

THEN  USED  TO  FIT  KITH  THE  E IP  £3 Id  ENT  AL  DATA  TO  IDENTIFY 

J 

J 

THE  CATHODIC  AND  ANODIC  REACTION  RATE  CONSTANTS  OF  THE 

c 

c* 

SURFACE  REACTION. 

<w. 

c 

V. 

c 

-  NO d EH CL A. USE  - 

- 

N- 

N.(=:.UHBEG  o F  PoINTS  IN  X  DIRECTION. 

,t 

N=NUM3ER  Or'  INTERIOR  COLLOCATION  POINTS. 

c 

r* 

NT= NUMBER  OF  TOTAL  COLLOCATION  POINTS  INCLUDING  TNG 

c 

c 

JCUNDitaY  P Oi NTS. 

n 

v- 

1(1*1)  ,  Y  { 1  ,  L )  ,  .  .  .  Y  ( 1  ,  N)  =  CON  CENT  5  A  TIO  NS  OF  NITRATE  ION  AT 

Z 

r 

S  INTERIOR  COLLOCATION  POINTS. 

z 

v» 

Y  ( 1  ,  :<♦  1 )  ,Y  (  1 ,  &  *2)  ,  .  .  -Y  (1 , 2*N)  =CONCEN  TR  AT  10 NS  OF  H  YDa0.vY 

c 

c 

ION  AT  N  INTERIOR  COLLOCATION 

c 

u 

POINTS. 

V- 

c 

Cl  (1)  *  Cl  (2)  cl  (NX)  CONCENTRATION  OF  NITRATE  ION  AT  EACH 

c 

c 

POINT  IN  X  DIRECTION. 

c 

s» 

C2  (  1)  ,  .2  ( C)  ,  .  ,C  2  (NX)  CONCENTRATION  OF  HYDROXY  ION  AT  tAC.i 

c 

c 

POINT  IN  I  DIRECTION. 

c 

\  z 

X  S  ( 1 )  3  3U  R  F.iC  o  CO  NC  E  N  TR  AT  I CN  OF  NITRATE  UN. 

c 

C 

X3  (2  )  3  SU  R  ?  AC  -  CO  N"C  E  N  TP  AT  I  ON  or  HYDaoXY  ION. 

J 

!  c 

C 1 0=  BJ LX  CONCENTRATION  OF  NITRATE  ION. 

z 

t 

C2  3=  3li  LF.  C  JNC  ENTRATI  CN  OF  HYD30CY  UN. 

’w 

c 

D'1AX  =  0  IF  FUSION  THICKNESS,  Cd. 

c 

1 

IS 3  LOC  AI  ION  Jr  THE  SPLINE  POINT,  0  <  S3  <  1  . 

c 

1  c 

02  3  =  i.N  CR  Ed  EN  T  OF  THE  SPLINE  POINT,  EG. 

c 

II  MAX3  MAXUUG  .IME  TC  BE  INTEGRATED,  SEC. 

c 

JT  Id E=  Tide.  INCa E.’IENT,  SEC. 

w 

Z 

T  A  UMAX3  D  I  d  E  NG  Io  N  LESS  MAXIMUM  TIME  To  BE  INTEGRATED. 

,'3 

l 

4  rt  a*  a*  O  a  1  '1  w 

;; 

c 

»  >  .\  J  —  J  1 .1  ~  ,i  J  2  sj  »«  L  Z Z  S  I  I  *1  Z  I  !IC  U  Z  kj  Z  r 

■J 

z 

irty  iivj.i  3 ") Z 3  C i’  ."i  L  D 2  J .(  i  1 0 v* • 

J 

1 

G  o  T  Ass  r>ACT  i  o  .<  u?JE?  CF  NITRATE  ION. 

f  ^ 

c  o  i ii  j « *  i  «*  o  Jo  cl  u a w «.  r o j  l\ a i u  <«j .< oiu  a  t,  x ,\.i  j  jio 

J 

Z  *.  ii  w  *  1 G  it  • 

c 

- 

OO.i  1  =  H  ol  oR  D  G  -  -i  wO  US  REACTION  RATE  CUNGTAnT  IN  CAT.iJD.C 

Z  a  *t  wC 1 1  w  »•  • 

CD-CURRENT  OE.iSTTY,  Ad  P/C  d*  •  z  . 

-i 

- .  • 
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THE  DIEFJSIJ.i  COEFFICIENTS  ■>:  ISC  NITRATE  A.<D  HOWH 
IONS  ARE  1. <#022-5  A  N  D  5.  2oOE-5C22/SEJ,.R  E3PECTIV  EE  I. 

i*********************«*’ 


v. 

C 


c 

*c 


IMPLICIT  REAL'S  (A-H,  0-2) 

DiaENSUN  A  (20,10)  #a  { 2  □  ,2  o)  ,C  1  (oJ)  ,C2  (60)  ,C  ( 2  0  j  ,  I  (6  ,20)  ,  I  .'.AX  <2  O)  , 

1  Dir  1  (20)  .DIF2  (2  0)  ,DIr  3  (2  0)  ,  ROOT  (20)  ,  VECT  (20)  ,  ERROR  (20)  , 

2  S  A/2(24 , 20)  ,  ?N  (400)  ,  XS  (2)  ,  XX{  1)  ,  WA  l  J)  ,3A5  (1) 

EXTERNAL  A  J  Xe 

CC.1.1GN/L  1/C  J2E  1 ,  C  o  ZF5 
CO/.;iGN/L2/OOEF2 
CCil.lUN/LJ/N  ,  N  r 

co;lio.'«/z.4/a  ,  £ 

COiiaON/L5/\S 
COilEON  /Lb/COr.P  4 


C  U  .El  ON/  7  /  .i  S . 


no: 


C  0 .1 A  J  N  /  L  d  /  C  O  E  E  3  ,  C  O  N  1  ,  C  G  N  2  ,  O  A  2  A  ,  3  S  £  A 

cg;«cn/l9/x 


I N  P  J  4.  DATA  -■ 


READ  (5 , 3  10)  N D,  N,  NO,  N1 , AL,  3 E, EPS, TIE A  X,  IC, DOS 
AEnJ  (5,320)  C1S,C23,D/iAX,DTj..'lE,Nv,ES 
READ  (5  ,  S  J  0)  GArfA,  B  ETA,  CCN2,  CO. i  1 
B  EAD (5,340)  j  S 10 , N  EE ,E  P 

- CALCULATE  SO. 'IE  CONSTANTS  TO  oZ  USED  IN  ONE  PSDSRA/. - 

C0EF1  =  ( 1. 0  J2 0- 5)  / 5.  2  b  D  -  5 
CO.i ;  2  =  C,.  9/  C  1  £ 

COE ?3  =  ( 1 .  J02D-5)  *C1 S/D 1A X/ES 

COE  F4  =  2  .  *  (  1. 902D-5)  *C1  3*9o43  7./0/.AX/ES 

COEr5=1./4,ci/ 

D  a  =  D  iIA  X/  (Ni-1) 

DT AN  =  (  a.  2  o  J  -  5)  'DTI  XE/0  3  A X/D  .1 A  X 

T  A  J  .1  A  X  =  ( 5 .  E  o  j-  5)  * T  IX  A X  /  D.l  A  X/  D .1 A  X 

SDTA J=uTA  U 

DETA= 1 .  /  (NX-1) 

nt=n  ♦no*:;  i 

NE3=/*N 

- SET  INITIAL  DI  R  EN  S  I  ON  4.  E  SS  TI.1Z  =  0. - -  — 

I A J=  0. 

- 0SZI..E  7  ii  L  I'AaIITJO  A.iJ  .11  N  l:U  0  D  111  ENS  I J  N  LES 

.  i  4  I  N  C  4  4fc.  N  4  4  0  .4  M  J  .4  W  0  i  N  4  .1  .1  4  .1  4  4  J  tl  .t  4  .U  .1 

SUE  joU  I.  N  w  D  . r  j  ii 5  —  —  —  —  —  — 

OTAJ  .11  =  .  0  0  1  *  L  T  A  U 
DTrtU.IA  =2.0*  OTA  J 
«  ?  I T  E  ( u  ,  j  j  0  i 


1 


n  n  r,  r,  n  n  n  n  n  n  r,  n  Ci  n  n  r.  n  n  n 


140 


- E  V  A  LJ  AT  E  THE  COLLOCATION  POINTS  (BOOTS)  OF  THE 

JACOol  POLYNOMIAL  OF  OLDER  N,  AS  N  ELL  AS  T ci 
FI  a ST  AND  SECOND  OSH I V AT  IV  35  OF  THE  POLYNOMIAL 
hT  TaE  HOOTS.  - 


CALL  J C 0 0 1 l N D , N , NO ,  N 1 ,  A  L , 3  S ,  D I F 1 , 0 1 F2 , DIF3, HOOT) 

WHITE  (6 ,350  j  (HOOT  (J)  ,  J=  1  ,  NT) 

- CALCULATE  THE  D ISC A  ET ISA II ON  COEFFICIENT  MATH  1 1  A - 

10  =  1 

DO  20  1=1, NT 

CALL  DFO?  d  IND,  N  ,  NO  ,N  1  ,1 ,  ID,  OIF  1 ,  DIF  2,  3173,3001 ,  '/EC  I; 

DO  10  J=  1, NT 
10  A  (I,  J)  =  /ECT  (J) 

20  CONTINUE 

- CALCULATE  THE  DISCRETISATION  COEFFICIENT  idli til  0 - 

ID=2 

DO  4  0  I  =  1  ,  NT 

CALL  D  F  0  ?  it  ( M  0  ,  S  ,  NO  ,  N  1  ,I,I0,DI?1  ,  DIF2  ,  DIF3  ,  HOOT  ,  V  EC.  ) 

DO  30  J  =  1 ,  NT 
30  3 (I, Jj  =  VSCT  (J) 

40  CONTINUE 

- X N p J x  THE  VALUES  OF  ARGUMENTS  TO  5E  USED  III  THE 

INTEGRATION  SUBROUTINE  DIF3U3  - 

Mr  =  1 

JSTA  HI =0 
MAXD  S&  =  7 
DO  4  5  1=1, N  EH 
IMA  X  ( I)  =1.0 
45  CONTINUE 
1  =  0 

-  INPUT  THE  INITIAL  CONDITIONS - 

DO  5  0  J  = 1 ,  N 
Y  (1  ,J)  =  1. 

I  ( 1  ,  J  *N)  =CO£F2 
50  CONTINUE 
XS(1)  =  1. 

X<j  U  J  * .0 . F  ** 

WRIISfo,  3*0) 

SO  1 1  ME  =  I A  U  *  D  .1  A  a  *  D  M  A  X  /  (5.  2  60-3) 

rtiiD  ?  R  I  N  .  i  H  r*  C  J  H  3  c*  N  .  D  r.  N  o  I T  i  AT  k  ;i  I  j  I  M  r, 
CONCENT  5  A  TUNS  OF  NITRATE  AND  .i'iZ  S  0  X  i  lo  .,5 
COLLOCATION  POINTS - 

CALL  C  J  3  S  N  T  \,C  D ,  A,  Y  ) 


™  Caii.^.aLAAn 

Pa  IN.  THE 
A  a.  .HE  NT 
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c 

c 

c 


X 


c 

c 

c 

c 


*RITE(o,3oO)  TI H E, CD,  X 3(1)  ,  X  a  ( 2 )  ,  K  r  L  A  G  ,  I T  S  R  N  ,  a 
Wair£(S,70J)  (  Y  (1,1)  ,!  =  !,.<  EE) 

-  iVnxJATS  AND  PRINT  THE  CO H C ZN  T  RATI C NS  Or 

nil  H  AT  Z  AMD  HY  D  30  X  'l  IONS  A  2  EACH  POINT 
IN  A  DIRECTION  AT  THIS  TIHE - 

CALI  CON DEN  (C 1,C2, NX  ,N  D  ,D  ET  A,  ROOT,  DIF  1  , Y) 

WRITE  (6, 230;  (Cl  ( J)  ,0=  1,  NX) 

WRITS  (6,2 -jD)  (C2  ( J)  ,J  =  1 ,  NX) 

- C1ECX  IS  THE  INTZSA  AilON  TI.1S  AT  THIS  .1 0  A  ENT  a  £  A  C  n  S3 

Itix  H  A  X  HI  U  '1  TINS  TO  32  INTEGRATED  — — - 

IF  (TAU.GE.  TAUHAX)  GO  TO  1  >0 
85  !=!♦  1 

IF  (I .  E  J.  1 )  Gu  TO  9  0 


x 

C 

C 


IEST  IF  THE  CONCENTRATION  G?  NITRATE  ION  AT  THE  NTH 
COLLOCATION  POINT,  1 2,  TH2  LAST  ONE  SEEDS  E  THE  SPLINE 
POINT  IS  WITHIN  0.CC01  Ox HEN SION  LESS  CONCENTRATION 
JNif  OS  THE  BOUNDARY  CONDITION.  - 


C 

C 

c 


IF  ( (  1. -Y  (  1,  J)  )  .  LT.  0.  0001)  NO  TO  30 

-  IF  TaS  TEST  SAILS,  THE  SPLINE  POINT  *3  INCREASED, 

HOVED  FURTHER  INTO  THE  SOLJTION. - 

CALL  CHANGE {CS, DZ5, CO  Z  S3  ,  C  u  Z  S  C ,  C  J  Z  S  5  , FACTOR, D  T  A  J ,  D  T  A  J  HI , 
1  J  START) 


D  I A  J  H  A  , 


EVALUATE  THE  CONCENTRATIONS  OF  NITSATs.  AND  HYDROxY 
1 0 N  S  AT  THc  CO  L  LOC  ATI-/  N  PO I  N  T S  AT  T  H  x  N  xW  COO R  D I N a  T  E 
OUx  TO  THE  INCREASE  0?  THE  SPLINE  POINT.  - 


C 

CALL  EXPAND  l Y, BOOT, DIF  1, FACTOR, ND) 

C 

c  - integrate  equations  { 5 o d  and  ;3i)  using  the  concents atix 

C  Ax  Tna  wOL  xG  CATxON  ?  0  x  N  .  S  AT  THIS  T I H  x  TO  0  &  T  a  x  »*  THE 

C  CONCENTRATIONS  AT  THE  SUCCESSIVE  II 'E  - 


9u 


C.ALx  D  IS  S  U  i.  { U  ER  ,  T A  U,  Y,  S  A  VE  ,  DTA  J  ,  DT  AUU  I  ,  0 T A  U  HA 
x  3  3  0  5  ,  X  F  LA  G  ,  G  S  T  A  3  T  ,  .1 A  X  0  E  3  ,  ?  N ) 

IS  (XSLAG.  LT.  1)  SC  TO  190 
x  x  rtU  1x=.  JO  1  *  0  i  A  U 
JT  A  U  H  A  =  2  -  0  *  0 1 A  U 
xF(x.Gx.x)  jO  TO  9  3 

IF  U  (1,1)  .0  1.  1  ..  OP  .  Y  (1,2)  .or.  1 - .3.  Y  (  1,3)  .01. 


?  j  ,  1 1  r  /  .  H  a  a  , 


.03. ili, x)  *  ,  x  .  i  . 


n  n  n  o  r  i 


1 


1>2 


.  OR.  OAoS  (  T  ( 1 , 5)  -  1.  )  .  G  r.  0-  0304)  JU  TO  90 
GO  TO  93 
9  3  :n=i/ic 
M2=M  1*IC 

IF  (I.  E  .12.  OS.  TAU.  GE.TAUKAX)  go  TO  95 
GO  TO  db 

-  CALCULATE  THE  SURFACE  CONCENT RATIONS  - 

95  XX (1) =Y (1,  1) 

ITSa  N=  50  0 

CALL  ZSXSTM  ( A  U X2 , S P, NS  13,  NEC,  XX  ,1 1  ESN  ,  •  A  ,  B  A  LER) 

XS  (1)  3  XX  {  1 ) 

XS  (2)  =  -i.  *CJEF1*  (A  (  1,  1)  *XX  (1)  *A(  1,  NT)  ♦1.0)  /A  (  1  ,  1}  -  A  (  1  ,  NT)  /A  (1  ,1 ) 
1  *CJc F2 

DO  150  J=1,N 

XS  (2)  =  XS  { 2)  -  3.  *CCEF1*A(1,J*1)  *Y  (  1 ,  J)  /  A  (1,  1)  -  A  (  1 ,  J  +  1 ) 

1  «*(  1,  J*N)/A{1,  1) 

150  CONTINUE 
GO  TO  60 

190  WRIT  2(b,  240)  DX,S0TAU,  N  X,  DT  AU  ,  DTJ.  ME,  DM  AX  ,  AL,  o  5,  ZS , EPS, X FLAG, 

1  E? ,  NS IG 

NnirZ(0,2‘O)  SETA,  GAMA, CONI  ,CJN2 

-  FORMATS  FOR  INPUT  AND  OUTPUT  STATEMENTS - 

2  30  FORMAI (/,26X, *N03- •, 1X,9D  1  1.4/30X,  901  1.4/30X,  90  1  1.4/30X,  9  0  1  1. 4 

1  /3 OX ,901  1. 4/30 X,  90  1 1.4/3 OX,  901  1.  4/30X, *0  1 1.4/30X,  *0  11.4 

2  /J OX, 901 1. 4/3 OX, 90 1 1.4/3  OX,  9011. 4/3  OX,  *0 1  1.4/30  a, 90  11.4) 
240  FORMAT  (//.  5X,  '  X  INTERVAL3'  ,012.4,32X,  '  INITIAL  TAU  STEP  Sj.ZE=' 

1  ,  01  2.  4/5X  ,  '  NO  OF  POINTS  IN  X  DIP*  <,  13, 2'Jji, 

2  'LaST  TAU  STEP  SIZE3  '  ,  0 1 2 . 4/5  X ,  '  RE  AL  TIME  INTERVAL3  ' 

3  ,D10.4,24X,  'DIFUSION  LENGTH3  '  , 0 1 0. 4//5 X , ' AL3  ' , F 7 . 2 , 3 OX , 

4  '  u  E  =  ' ,F7. 2,25X,'ZS3  ' , F7 . 4//5X, ' EPS3  ',D10.2,27X, 

5  '  K  F  L.kG3  '  ,  I3//5X,  '  E?  3  '  ,  01  0 . 2 , 27  X,  •  NS  I G*  *  ,  I  j) 

245  FORMAT  (//,5X,  'REACTION  ORDER  G?  N03-  ION  =',F7.2,25X, 

1  *  a  E  ACT  ION  ORDER  OF  OH-  ION  =  '  ,  F7. 2//5  X ,  '  1  ST  CONSTANT  3  ' 

2  ,j12.4,25X,  '2ND  CONSTANT  =',012. 4) 

2  50  FORMAT (/, 25  a, 'OH- ',21,9011. 4/3  OX, 9011. 4/3CX,  9  0  1  1 . 4/30X, *0  11.4 

1  /JOa, *  Oil.  4/3  OX,  90  11.4/3  OX,  901  1 . 4/30  X,  9  0  1  1. 4/  33  X  ,  9  0  1  1 . 

2  /3Ua,901 1. 4/3  OX, 901 1.4/3  OX, 901 1 .4/3  OX, *0 1 1.4/30X, *011.4) 
310  FORMAT  (413, 2F5-  1 , 201  3.  2  ,1 5 ,  Fd.  4) 

3  20  FORM  AT  (4012. 5, 17,  F7.  4) 

3  30  FORMAT  UF7.r,  2012.  4) 

3  4  3  FORM  AT (215, 01 2. 4) 

3  50  FORMAT  ( 1  5  X  ,  1  0  F  1  0  .  6  ) 

3  60  FORMAT  (//,  0  1 „ .  4,  1  X  ,0  1 2 . 4 , 2  a.  ,2  0  j  0.  1  J  1 5  ,  F2  0 . 7 ) 

37  0  FoR  MAI  (/,2a,oF2U.  6/2  X,  oFJO.  oi 
3  90  FORMAT  (2X, 'CJLLCCATICN  POINTS:  ') 

3  90  FORM  AT  l//,  3  a,  'TI9.E  ',  12  X  ,' CURRENT'  ,4bX,  '  ION  CO  NCE  NT  LATT«.i  '  ) 

7  00  FORMAT  (3  OX,  1  0  F  1  J  .  5/3  Ja,  1 JF  1  0.  5) 

STOP 

END 
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SUBROUTINE  JC03I  (N  D,  N,  30,  N  i  ,AL, 32,  DI FI  ,  DIE  2,  JIF3  ,  ROOT) 

I  M  PL  I C 1 1  hh  X'l*  6  (  A  -  H  ,  0-  2  ) 

DIMENSION  D  iF  1  (ND)  ,DIF2  (ND)  ,DIF  3  ^  1«0)  ,  LOOT  (ND) 

5U3R0UTIN E  EVALUATES  THE  COLLOCATION  POINTS  (SOOTS)  CF  THE 
J A CO 31  POLY  NOMI AL  OF  ORDER  X  AS  N ELL  AS  THE  FIRST,  SECOND 
AND  THIRD  DERIVATIVES  OF  THE  POi.Yaj.UAL  AT  THE  ROOTS. 

- IN  PUT  PARAMETERS - 

INTEGER  i'i  D=  T  H  E  DIMENSION  OF  VECTORS  Di  F 1  ,  3IF2  ,  D I  ?  3  ,  ROOT. 
INTEGER  N  =  iHE  DESS  EE  OF  THE  JACOBI  POLYNOMIAL,  I.E.,  THE 
NUMBER  OF  INTERIOR  COLLOCATION  POINTS. 

INTEGER  NO* DECIDES  WHETHER  X=0  IS  INCLUDED  AS  A  COLLOCATION 
POIwT.  NO  MUST  3  E  SET  E  2  J  A  L  TO  1  (INCLUDING  1  =  0) 
OR  0  (EXCLUDING  THIS  POINT). 

INTEGER  N 1 =  A S  FOR  NO,  BUT  FOR  THE  POINT  1=1. 

REAL  An,  3  2=  Til  Z  PARAMETERS  OF  THE  J  AC  031  POLYNOMIAL. 

-  OUTPUT  PARAMETERS  - 

REAL  ARRAY  EDO T=ONE- DI ME N SI  ON AL  VECTOR  CONTAINING  ON  Ell. 

THE  JN-N0  +  N1  DEEDS  OF  THE  POLYNOMIAL  USED 
IN  THE  COLLOCATION  ROUTINE. 

REAL  ARRAY  =  T  HR  E  E  ON  E-DIM ENSIGN AL  VECTORS  CONTAINING 

DIF1.3IF2, DIF j  ON  EXIT  THE  FIRST,  SECOND,  AND  THIRD 

DERIVATIVES  OF  THE  POLYNOMIAL  AT  THE  EEROS. 


;  ********  *  *  *  *  **£ 


C  - FIRST  -VALUATION  0?  COEFFICIENTS  IN  RECURSION  FORMULAS  - 

C  - RECuAjIQN  COEFFICIENTS  ARE  STORED  IN  DIF1  AND  DIF2  - 

C 

A  9=  A  L*  3  E 
AD=BE-AL 
A?=3  E*  AL 

D1F1  ( 1)  =  (AD/  (A3*2)  *  1 )  /2 
DIF 2  (T) =0 

IF  (N  ,LT.  2)  GO  TO  15 
DO  12  1=2, N 
-  1=1-1 
-=A3  *2*-1 

DIF  1  (I)  =  (AS  *A  D/2/  (  Z*  2)  *  1)  /2 
IF  (I .  NE.  2)  GO  TO  1  0 
DIF2 (I)  =  i»j*rt?*21) /Z/Z / (I  *■  1 ) 
ju  TO  12 
2  =  2*2 

1=21  *  ( A3*  -  1  ) 


10 


c 

c 


lc-ii 


t=i*  u *  +  'L) 

Oj.72  (I)  s£/„/ii-1) 

12  CO 21  IN J  E 

- - -  30^,1  DOT  27 I  NA  TI  0  N  3?  ME*  ION  MZTHOu  Wir.’i  3  'J3  ?F.  £3  31  J  2 

OF  PLaYIOUSLY  DET  Z3  7.1  3F  J  aU0T5  - 

15  X=0 

Du  35  1=1,2 
20  XD=0 
XN=1 
:<D1=0 
X  S  1  =  0 

t/O  2  x  J  =  1  ,  .i 

X?=  (0171  (J;  -a)  *XN-CIF2  (J)  *.c^ 
xpi=  (difi  (j)  -X)  *xin-oiF2(j)  *xdi-xn 

XD=X2 

xd  i=x:m 
X.\  =  Xi> 

22  x:n=x?i 

2C=  1 

z*x:;/xa 1 

IF  (I .  E*.  1)30  ro  3  0 
DO  2  5  J  =  2 , 1 

25  uC-ZC-Z/  (X-SOOT  (J-  1)  ) 

30  L=0/2C 
A  =  X-  - 

IF(DABS(-)  .01.  1.D-C9JGC  TO  20 
BOOT  (I) =  X 
X  =  X».00J  1 

35  continue 

- ADj  2  V  3  2 1  u  A  L  COLLOCATION  POINTS  AT  X=J  03  1=1 - 

NT  =  N  +  NO  4  ;t  1 

.r  (NO.iiiy.1/}  \jk)  ro  4  5 
DO  40  1=1,3 

J  =  .l  4  1  -I 

uo  sour  (J  4  1  )  =ii0u7  (J) 
nOOT  ( 1)  =0 

45  IF  (N  1.  2*.  1  )  NOOT  (NT  )  =  1 

- NO.»  27  A  L  J  AT  E  DEtiIVATI/Z3  OF  POL  I  30.11  AL - 

Du  5  0  I = 1 , 3  ? 

(I; 

DIFI  (1/  =  1 
0IF2  ( 1)  =J 
JIF3  ( X ; 

DO  o  u  J  =  1  ,  ,i  ! 

*  -  (j.  —  •  ij  j.  *  0  3  0 

i~x~ aod  r  ( j) 

DIF3  (Ij  =  F  *  0  1 1-  J  (I)  4  3=0172(1, 

DII2  (I,  3  I  *  J  IF  2  (I)  «■  2*  DIF  1  (I) 


t 


( i  n  <  i 


14 


3  IF  1  (I)  =  x  *D  IF  1(1) 
50  CONTINUE 

fztu  a:i 
end 


s JuaouriwZ  jf opr  (n d,  n,  no,  n  i  , i,  i o,  qifi  ,  oifi,  dif 3 , aooi,  v - c r ) 
1.1  PL.  10  IT  RSAL*d  (A-H,0-  3) 

3 1. IE  NS  ION  3IF1  (NO)  ,  DIF2  (NO)  ,  DIFS  (3  3)  ,3uOT  (HD)  ,  VSCT  (33) 


*******»»»*■ 


C  SOjkOUriilS  EVALUATES  DISCRZTIZA  "ION  COEFFICIENT  MAIJICZS  c 

<:  AN  3  GA  USE  j.  AN  QUADRATURE  WEIGHTS,  NC.1  ALT  ZED  TO  S  J  .1  1-  C 

c  c 

C  - INPUT  PAR  A.'l  2TE  S  S -  C 

c  c 

C  INTEGER  N,  JO,  N1=AS  IN  SUBROUTINE  CCC31.  C 

C  INTEGER  I  =  THc  INDEX  CF  THE  NODS  FOj  WHICH  THE  WEIGHTS  ARE  v. 

TO  ot  CALCULATED.  C 

INTEGER  ID=1N 3ICAT0R.  13=1  GIVES  THE  WEIGHTS  FOR  DY/3X,  C 

13  =  2  GIVES  THE  WEIGHTS  FOR  D2Y/JX2,  AND  13=  J  v. 

GIVES  THE  GAUSSIAN  WEIGHTS.  THE  VALUE  OF  I  IS  C 

0  la  RELEVANT  IN  THE  LAST  CASE.  C 

C  REAL  ASSAY  ROOT,  C 

C  JIF1,3IF2,  JlrJ  =TH  E  0 NS  -  01  HE  NS 1 0 N AL  VECTORS  CONFUTED  IN  0 

C  SUBROUTINE  JC03I.  C 

C  C 

C  -  OUTPUT  PARAMETERS  -  C 


C  REAL  A  a  3  AT  VECT=THS  COMPUTED  VECTOR  OF  « EIGHTS.  C 

C  0 

NT=N +N0+N 1 

IF  (I  D.  s;.  J )  Go  TO  1  0 
30  20  J  =  1  , J  . 

IF  ( J  •  N E.  I)  G 0  TO  21 
IF  (I  3.  nZ.  1 )  Gu  TO  5 
VECT  (I)  =Di?2  (  I)  /DIF1  (I)  /2 
Gu  TO  20 

5  VSCT  (I)  =JxFU(I)/DIF1  (I) /3 
GO  . o  20 


21 

Y  = SO  OT  ( 1 ) 

-LOOT  (0) 

VECT  (J)  =3 

IF  1  (I)  /3IF1  (J)  /Y 

xr  .13-  x * . 

2)  VZCT(J)  =  V  EC  T  (J)  *  (3IF2  (I) /JLF^ 

r  > 

CwNliNUE 

jO  4  3  3  0 

10 

l  =  0 

O'.,  2  5  J  =  1 

,i  r 

::  =  ROui  ij) 

A  X  =  X  *  (  1  -  a  ) 

1F{N0.S«.  J)  A  *  =  \  X  /  X  /  X 

(.1  L  4  ,  .  3)  «i  X  =  A  \ /  (  1  -  X )  /  ( 1  -  a  ) 


1L6 


V EOT  (J)  =AX/JIF  1  { J)  **2 
25  l  =  Y  ♦  /  E  G  T  [  J  ) 

CO  o  0  J  =  1 ,  N  I 
60  VECI  (J)  =  Vr.CT  (J) /Y 
50  HETURN 
END 


s UDRoariNZ  int bp  (s c, nt, x, soot, di?i , xintp) 

Iii?L  10  IT  S2aL*8  (A-H,0-Z) 

DIMENSION  SOoT  (  N  D )  ,  DI  F  1  (NO)  ,XINT?(ND) 


-■******, 


««***«*«* *» , 


c 

C  503 SOUTINE  EVALUATES  THE  LAGS  AN 01  AN  I  NT 


*******<1***********»****(J 

G 

ESP'JLATION  GgEFF IGIINT5 


C  -  INPUT  PAS AUSTESS  C 

C  G 

C  INTEGER  N  I  =  I  tt  2  TOTAL  N’KIjES  OF  COLLOCATION  POINTS  (  =  N+N0  +  ..1)  C 

G  FJa  WHICH  THE  VALUE  0?  THE  32P2NDZNT  V  A  SI A  RLE  Y  G 

C  IS  KNOWN.  C 

G  HEAL  X=TtiZ  A3SCISSA  X  WHERE  Y  (X)  IG  DESIRED.  G 

g  PEAL  ASS  AY  C 

G  SOOT, D I? 1  =COLLOCATICN  POINTS  AND  DERIVATIVES  CF  NODE  C 

C  iJu'i  N  0  P.  I  A  L ,  DERIVED  IN  GU5RO  U  TINE  JC'j3I.  C 

G  C 

C  -  OUTPUT  PiaAHETEHS  C 

•G  C 

G  SEAL  ASSAY  XI  NT  ?=T  H  E  VECTOR  OF  INTERPOLATION  WEIGHTS.  C 

G  G 


*******  i 


:«#*»♦*  *G 


?CL=  1 

DO  5  I  =  1  ,  N  T 
Y= X- ROOT  (1) 

XINTP  (I)  =0 

IF  lY.E^.  0.  00)  XINTP  (I)  =  1 
5  ?GL=  POL  *  Y 

IF  (POL.t*.  0.00)  GO  TO  10 
DO  o  1=1, NT 

o  XINTP  (I)  =POL/DIF  1  (I)/  (X-30oI  1 1)  ) 
10  SETUHN 
END 


jJJnuUi'u 

1 

T.-.p-icir 
j  The  ns  ion 
3 

0  x  X  g  4*  S  I J  N 


E  OikjJo  (*i|  »  ,  ’)  A  /  . ,  d  ,  >i>.  x  N  ,  *1  ■!  .C  ,  up  j  f  ,'i  r  ,  i 

JS TA  S  T ,  U  A  X 0 E  A  ,  P  N ) 

9  u  liL  ^  1  [  A  *  il,  ^  “  w  ) 

L  (0,20)  ,  Y.l  AX  (  20  )  ,  SAVZ(2-*,tJ)  ,  ERROR  i  E  J  ) 
A  (o  i  ,P'£PTS7(7,2,j) 

<.  UO)  ,?JL(  20 ,20)  ,X  (2  0) 


-lii 


, 


O  J  O  T 


,  ?«  k-t  0  J  j  , 


«T  Li  rt  VJ  f 


*  *  **  ******4rft«*«*att*v««*a*««*;j 


o  n 
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\ 


\ 

i 


i 


c 


S  U  xl  R  0  J  Cl  HZ  xNlZGRAT 
ORDER  EQUATIONS  AS 
ORE  STEP  0?  LENGTH 
DECREASED  .iliulN  TH 
AS  LaaGE  A  STEP  AS 
STEP  ERROR  '#  a  IC  !1  IS 

each  component  of  t 


ES  A  SET  0?  .<  ORDINARY  DIFFERENTIAL  FIRST 
DESCai  JED  i:i  E  *  U  A 1 1 0  N  S  150)  AND  (51)  OVEc 
a  at  each  ca^l.  a  hay  he  increased  or 

Z  A  .Ml  OS  a  M I  N  TO  a  H  A  a  IN  ORDER  TO  ACHIEVE 
POSSIBLE  * SIZE  NOT  COMMITTING  A  SINGES 
LARGER  THAN  EPS  IN  THE  L~ 2  NOSH,  NHxRE 
HE  ERROR  IS  DIVIDED  BY  THE  CO . IRONS  l.TS  of 


Y  MAX. 


THE  PROGRAM  CALLS  FOSS  SUBROUTINES  NAHSD 
DIFFUN (T,Y,DY) 

H  ATiN  V  (P  N,  N  ,  M  ,  J) 

SOLVE  (  N  H  ,  t* 'J  L ,  3 ,  X) 

PEDERV (T , Y , ?N , M ) 

THE  FIRST, DIFFUN, EVALUATES  THE  DERIVATIVES  OF  CHE  DEPENDiNI 
V  ARIAS  DCS  STORED  IN  Y(1,I)  FOR  1=1  TO  N,  AND  STORES  T  ;iE 
DERIVATIVES  IN  THE  ARRAY  CY.  THE  SECOND  IS  CALLED  ONLY  If 
Hr  IS  SET  TO  1  CS  2  FOR  STIFF  METHODS.  IT  ?ESF OR  IS  A  rlaJT 
STAGE  OF  GAUSSIAN  ELIMINATION  ROUTINE  OF  TiiE  N  3Y  N  MATRIX 
STORED  IN  THE  ARRAY  PM  (M , M)  .  SOLVE  PERFORMS  THE  jACK 
SIJSSTITUilON  PROCESS  0?  THE  GAUSSIAN  ELIMINATION  RoUTIN-. 

C  u  ^  o  ‘O  ji  j  u  2  u  J  0  il  u  Y  xr  Hr  To  1,  A. ID  ^  k  x  » l  £>  P  a  n  .  x  a  .« 
DERIVATIVES  Of  THE  DIFFERENTIAE  EQUATIONS  PSOV I.DED  DY  DIffJN. 

- PARAMETERS  — - 


N=  THE  NUMBER  OF  FIRST  OR  D  cS  OIFFERsN  UAL  EQUATIONS. 
1=1115  INDEPENDENT  VARIA3L2. 


c 

Y— 

A  d 

BY 

N  ARRAY  CONTAIN 

ING  THE  DEPENDENT  V 

ARIA3LES  AND 

c 

W 

THEI 

h  1>Z 

ALSO  DERIVATIVE 

S.  THE  DEPENDENT  VARIABLES  Ahc. 

c 

z 

S  *  o  £1 

ED  I 

N  Y  (1  ,1)  ,  1=1  T 

C  N. 

V- 

» 

w>  A 

V  E=A 

DLJ 

CX  Of  AT  L  £A  ST 

12  *N  FLOATING  POINT 

^GwAirCIi^  i  j  w  j 

c 

c 

b 

Y  TH 

E  SUBROUTINE. 

c 

H= 

THE 

3TZ  ? 

Si  Gw  TO  3  E  ATT 

EM  PIED  ON  THE  NEXT 

STEP. 

w 

HM 

IN=I 

HE  MINIMUM  STEP  SIE 

E  THAT  MILL  BE  USED 

FOR  THE 

0. 

1 

»i  r  z  j 

RATION. 

c 

:ih 

a  a—  r 

he  a 

A  a I M  U  H  SIZE  TO 

MHICH  THE  STEP  MILL 

3  w 

c 

E?S  =  Til 

Z  irt.T  Jit  *  Z  Z  r  CCNSTA 

NT- 

c 

J 

=  THE 

asr 

HOD  INDICATOR. 

Vw 

0 

AN  ADAMS  PREDICTOR  CORRECTOR  IS 

USED. 

c 

Z 

1 

A  MULTI-STEP 

METHOD  SUITABLE  FOj 

STIFF  S  Y  S  T  E  M  S 

Z 

IS  USED. 

c 

Z 

Z 

I  H E  SAME  AS  C 

ASE  1,  EXCEPT  THAT 

T  H I  S  S  U  3  R  o  U  T  x  N  _ 

z 

■* 

COMPUTES  THE 

PARTIAL  DERIVATIVES 

J  :  xi  KUaAw.u 

z 

£ 

DIFFERENCING 

OF  THE  DERIVATIVES. 

c 

z 

V  •/ 

A  X  =  A 

.1  A  a  \  i  ur  a  L  CC  AT  I 

0 N  S  Mu I CM  CONTAINS 

4.  .  t  X>  .!  .i  l  ll  j  M  V>  fc 

z 

j 

5 

ACH 

Y  SEEN  SO  FA?. 

z 

EZ 

ROR= 

AN  A 

R R  A i  0?  N  ELEME 

.<?  3 

X  F.w  X  40xlx 

c 

u  A  Z 

wPEt-  ERROR  IN  Z 

ACH  COM  PON  SN  I. 

c 

v  ^ 

M  .1  xl“ 

A  CD 

m p - e r  ion  code. 

^  xi  £,  A  .J  ^  Miiw  b  »  ^  o  JOwCMjbr  >  o  • 

z 

J  5 

r  4  a  r 

=AN 

l<?  r;  d  i  ca  tg  r 

« 

z 

z 

. . 

.<  D  ER 

=  ruz 

xjAaIMJ'I  d  e?  :  V  1 

::v z  :.uz  3.UJZJ  3 z 

V  _j  -  J  • 

V. 

c 

PM 

-  A  d 

uO  w«\ 

OF  rt  T  LEAST  N *  *2  r LOATING  pu I N T  LOCAT.ONS. 

z 

1 


************************(2 
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****************  *  *  **  ***** 


DATA  PEHiJi  /2.J, 4.5,7.  33j,  10.42,  13.7,17.15,1.0, 

1  2.0,  U.  3,24.0,37.39,  :>3.  33,  70.0  3,  o  7.  37, 

2  3.  0, o. 0, 9.  1  o7,  12.5 , 1 5. 3d  ,  1. 0, 1 . 0, 

3  1  2.  d,/4. 0,3  7.  d9,  53-53,  70.  08,37. 97,  1.0, 

4  1 .  ,  1.  , 0.5, 0.1667, 0.04133, 0.0  0d2t>7, 1.0, 

5  1.0, 1.0, 2. 0,1.0, .3157, .07407, .013 3/ 

DATA  A  (2)/-  1.0/ 

I H  ET  =  1 
Hr  LAG  =  1 

I?  ( JSTASI. LE. 0)  DC  TO  140 


1  oo 

DO 

1  10 

1=1,3 

DO  1 

10  J=  1  ,  K 

1  10 

GA  VE  ( J  ,  I)  = 

£ 

(0,1) 

HOLD  - 

a  3  2  « 

IF 

(a. 

EQ. HOLD) 

GO 

TO  130 

1  20 

£<  acj  ft 

=  11/ 3  OLD 

13 

ET  1 

=  1 

GO 

TO 

750 

1  30 

ilQOL  D 

=  y  w 

T  Oi.D  =  I 

k\  A  0  U  *1  —  i  .  0 

Ir  (JSTAdT. GI.O)  GC  TO  250 
GO  TO  170 

140  IF  (JSTA.-.T.  E*.-1)  GC  TO  1o0 
-«Q  =1 
:i3  =  n 

;n  =  n  -*  i  o 

N2  =  Ml  +  1 

;;  4  =  N  **  2 

:<5  =  n  1 

:i  o  =  ii  5  *  i 

CALL  DIFFJM  (T,Y,3A7E  (N  2 , 1 )  ) 
DO  150  I  =  1,0 

1  50  i(2,I)  =  EA7E(N1+I,1)  *.i 

:uiea  =  J 
K  =  2 
GO  TO  100 

160  IF  (  My.  E  iJ  vO  Lu)  JSTABT  =  1 
T  =  TOLD 
MO  =  N^Oi.0 


K  - 

♦  i 

GO 

TO 

120 

1  70 

IF 

(IF 

•  c*  0  • 

■j  0 

TO 

1  00 

1  F 

(  Ou 

-  g  r  • 

=> ) 

.j  G 

TO 

1  30 

GO 

TO 

(22  1 

■) 

/  «• 

22, 

223 

,224  , 

22  3,  22  o)  , * 

1  30 

L? 

( 

•  J  A  • 

n 

GC 

TO 

1  30 

AjO 

TO 

(2  11 

0  *• 

12, 

213 

,214, 

215,216,217) 

1  30 

Kr  L 

A  j 

=  “  *1 

r.  zz 

•J3:i 

2  1  1 

A(1)  = 

-  1. 

J 

GO 

:o 

2  30 

2  12 

a  r- 

)  = 

5  J 

J00 

00000 

149 


A  iJ|_ 

=  -  •) .  5  u  j  0  00  0  0  0 

JvJ  4 

u  3  3  J 

u 

13 

A  (1) 

=  -0.  -t  loo  6 06060006 6b 7 

A  l  3) 

=  -0. 750000000 

A  (‘0 

-  -j.  1  o  ob  6  D  o  6  6b  66  b  c7 

jO  T 

0  Z  30 

2  14 

A  (I) 

=  -0-375000000 

A  (3) 

=  -  0- 9  1  06 b  co  6  6 666  6  66  7 

A  (4) 

=  -0.  3333323333  33  3  33  3 

A  (5) 

=  -0.  0416666666  66  6  66  60-3  7 

30  10  230 

2  15 

A  ( 1. 

=  -0. 34oo  1  1  1  1  1 1  1  1  1  1 1  1 

A  iJ) 

=  -  1 .  0  4  1  66  6  o  5  66  66  7 

A  ( -) 

=  -0. 4jo  11  1  1  1  1 11 1  1  1 11 

A  p) 

=  -  o .  10416  6  6066666  66  7 

A  (o) 

=  -O-OOo  J333  3  33  J3  333J3  J 

jO  r 

0  2  3  ^ 

1 

1  6 

A  (  1) 

=  -  0. 32* 36  1  1  1  1 1  1  1  1  1  1  1 

A  13) 

=  -1.1  16066660060  66  7 

A  14) 

=  -J.o2o000G0C 

A  (5) 

=  -0. 177 )3 33 333333 33 33 

A  (o) 

=  - J. JziOOOOOCO 

=  0.  00  1  3  3  3  3  6  3  S  d  d  3  0  3  3  5  0 

Go  : 

0  230 

“> 

<- 

17 

All) 

=  -u.  3  1o53  13  3  12  16931  2 

A  l  3) 

=  -  1.  ZZ50J000C0 

A  (4) 

=  -0.  7  5 135135  1351 3516 

A  (5) 

=  -d. 255206333333  33333 

A  («) 

=  -0.0 436 1111111111111 

A  (7) 

=  -  0.  0  04  36  1  1  111  11  1  1  1  1  1  1 1  1 

A  (o) 

=  -0. 0  00  10341  26934  12o634 

oO  r  J  230 

2 

2  1 

A  tl) 

=  -I.OOdO'OOOOuO 

inO  *. 

0  230 

222 

A  (1) 

-  -0. 000606666660606 7 

A  13) 

=  -u.  3j J3J333333333.33 

r*  \  •** 
J  * 

O  23  ) 

■) 

23 

A  (  1) 

=  -d.  54545454  54  54  5-04 5 45  j 

A  (3) 

=  Ail) 

A  (4) 

=  -  0.  U0090303G9090  00  *1 

3o  r 

0  230 

T 

24 

A  (  1) 

=  -O.40OOOOOOG 

A  ( J ) 

=  -0. / OdOOOOOO 

A  1 4 ) 

-  -0. 2d0000000 

A  P) 

=  -u. J200000JC 

JO 

0  Z  J  j 

z5 

All) 

=  -0.  -*37  35 0 2 0  4  3  7  6  5  c  2 

A  (3) 

=  -  ).3Z)1o7jd22  11o73J 

A  14) 

=  -j.  3 lu216a731021 J4j 

A  (  5) 

=  -  0 .  J  5-»  7  44  5 2  55 -•  7  4 45  2  c 

— 

-3  •) 

-4  .  -t  1 6  3  5 J  36  4  36  j  5  J  4 

JO 

10  z  Ja 

1 

2  o 

A  ( 1 ) 

-  -j. -jo  1 6 32 5 53 Jc 1225 

A  (J| 

=  -u.  *  zG 0  3  4  5  2  Go  34  9  Z 0  0 
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A  (4)  =  -0 .  4  1  o  6 606  5  66 oo  o utj  o  7 

A  (5)  =  -0-  0  3*20476  1904761 5 

A  (fc>)  =  -0.  0  1  193476  1904  761  3 

A  (7)  =  - 0 .  0  00  566  3 9  34  240  Jo 2  3  2 
2  30  X  =  MQ  +  1 
I DuU B  =  K 
;ity?  =  (4  -  ;i?)/2 
EMQ2  =  .  5/FL GAT  (S3  ♦  1) 

BNOi  =  .  5/F  L J  A  T ( N  0  ♦  2) 

ZN*1  =  . 5/T  LOAT (UQ) 

PEPBTi  -  ZP3 

ZU2  =  ( ?  Z  3 1  a  i  ( !J  Q ,  M  TY  ? ,  2)  *  ?  E  23  d )  *  *  2 
Z  =  (PE3IST  (AC  ,«TYF,  1)  *  PE  P'Jii)  **2 
ZD2A  =  (2EoI3T  (AO,  .HYP,  J)  *2i?SH)  *»2 
IF  t2D«M.2Q.O)  GO  TC  730 
BN  Q  =  ZP  3*a  30  3/D  FLCA  I  (  N  ) 

2  40  I'AZVAL  =  111 

go  ro  (250  /O30)  ,  ihzt 
2  50  r  =  r  ti 

DO  2  bO  J  =  2 ,  X 
DO  260  J1  =  J,X 

J  2  =  X  —  J  1  *  J  -  1 
DO  2c  0  I  =  1,:,' 

2  60  i  ( J  2 , 1 )  =  Y  (J  2  ,1 )  *  Y(J2+1,I) 

Do  2  70  I  =  1,3 
270  2ac(0a(I)  =  0.0 

DO  430  L  =  1,3 

CALL  DIF FUN  (T , Y ,S A  VS ( 32 , 1 ) ) 

IF  (IW2VAL.  LT.  1)  GC  TO  350 
IF  (3F.L2.2)  GO  TO  310 
CALL  2EDZ32  (?,  Y,  FA  ,113) 
a  =  a  ( i )  *  d 
DO  230  I  =  1 ,  M 4 
280  ? »J  (I)  =22  (I)  *a 

290  00  300  I  =  1,3 

3  00  2  W  ( I  *  ( a  3  ♦  1 )  -  3  3  )  =  1.0  +  ?'.i  (I  *  (Hi*  1 ) -33) 

I  *  EV  AL  =  -  1 

CA  ul,  BAT  I  a  V  N  ,  3  3,  J1  ,  POL) 
i.F  (J1.GT.0)  GO  TO  3  50 
o O  TO  440 

310  DO  320  I  =  1,:i 

320  SAV^(^,I)  =  Y(1,I) 

DO  3  4  0  J  =  1 , a 

a  =  22  S*  3.1  AX  1  (EPS,  DA  35  (SAV£(4,J)  )  ) 


Y  (  1  , J)  =  Y  (1,  J)  ♦  ? 
J  =  A  (  1 )  *;i/n 
v.A  Lw  D  c  c  'J  ,1  (  .  ,  ( ,  j  A  / 
DO  3  j  0  I  =  1 ,  a 

z  ( a  6 ,  i )  ) 

3  30 

2  *  ( I  «■  (J  -  1 )  *  M  3)  = 

(3  A  7  £  ( :i  5  ♦  1 ,  1 )  -  D  A  V  Z  (A  1  *  I  ,  1 ) 

/  *■> 

3  40 

Y  ( I,  0  j  =  BA  7  2  (  5,  J) 

GO  To  2  .> o 

3  50 

it  iBF.az.O)  go  ro 

DO  3  o  0  I  =  1  , !( 

370 

360 

BA;B(*,I)  =  Y  ( 2 ,  I 

)  “  in  /  u  ( li  I  ♦  I  ,  1 )  *  :i 
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30 

ID  -4 

1  J 

3  70 

DO 

3  dO 

I  = 

1  ,  s 

3  80 

JA  7a 

{.45  +  1 

,  1) 

=  Y  (2,  1)  -  SAVE  {0  UI,  1)  * H 

call  s 

Da  7a 

(•< 

,?JL,5AVE  (N  5+1,1)  ,2) 

DO  399 

0=1  , 

.4 

3  99 

SAVE  (9 

,  0)  ** 

(0) 

4  10 

NT 

=  N 

DD  4  20  I  =  1,'l 

1(1,1)  =  i  (1.I)+A(1)*3A/S(9,I) 

7(2,1)  =  T  (2,1)  -  SAVZ{9,I) 

23*02(1)  =  ERROR  (I )  +  SAVE(9,I) 

1?  (DAJS  (SA  VE  (9,  I)  )  .  La.  (3:iD*Y.1A.C  (I)  )  )  t.  T  =  N  T-  1 
4  20  03. J  TIN  02 

IF  (OX.LE.O)  GO  TO  4*0 
430  CONTINOa 
440  T  =  TOLD 

If  (  (H.  La.  ( 1. 0000  1  ) )  .  A  ND.  (  (IWaVAL-.TXYP)  .  LX. -1)  )  uo  T 
IF  (  (Sc.  tw.  0)  .03.  (I«  EVAL.  :<£.  0)  )  3ACUM  =  a  ACU  A*  .  2  5D0 
IWaViiL  =  :i  c 
I3ET1  =  2 
30  TO  7‘jO 
^  60  KFLAG  =  -3 
4  70  DO  d  0  x  =  1  , 

DO  4  d  0  J  =  1  ,  K 
4  60  Y(J,xj  =  3A  VS  (J,I) 

F.  =  HOLD 
JQ  ~  NOOaD 
JSTART  =32 
RETURN 

4  90  D  =  0-  J 

DO  5  00  I  =  !,.i 

5  00  D  =  J  +  (a 3 FOR  ( I )  /  Y  M  A  X  ( I )  )  *  *  2 

I  *  a’/  AL  =  0 

IF  ID.  OT.  a)  00  TO  540 
I?  (rC.LT.J)  00  TO  =20 
DO  5  10  J  =  3,  K 
DO  5  10  1  ~  1 , N 

5  10  Y  { J  ,  i )  =  Y  ( J  ,  I }  +  A  (J)  *a&303  (I) 

5  20  XFLAO  =  +  1 
KNEW  =  II 

IF  (  xD  J  U  3 .  a  a.  1 )  30  TO  550 
IDOU3  =  ID0J3  -  1 
IF  (  IDOUo.  3  X.  1)  GO  TO  700 
00  5  30  I  =  1, 11 
530  3*  7E  {  1  0 , 1 )  =  Z3SCF.  (I) 

3  0  ID  /  0  J 

540  XFLA  j  =  i\:  U.lu  ~  2 

IF  (  a.  wi.  (JAIN *1.0 0001)  )  Go  TO  740 
T  =  TOLD 

xr  ( KFuAu*  u i«  *5)  jC  TO  720 
5  50  PP2  =  iD/a;  **2';i2*  1,  2 
Pa  3  =  1 .  a  *  x  0 

x  c  {  [  j  J  *  j  u  •  .1  n  a  D  a  E  )  mD  ft  .  (  Kc  ma^j*  yb  •  ~  1  j  j  'jo  TO  3  7  J 
D  =  0.  0 


0  4  oO 


1 


152 


DO  5  oO  1=  1  ,  N 

560  D  =  J  *•  (  (2350ft  (I)  -  SA  72  [1  J,I)  )  /  Yii  AX  (I)  )  **  2 
?SJ  =  (D/ZJ?)  **EN,}3*1.  4 
570  Pal  =  1. 

IF  (NQ.  L£.  1 )  GO  TO  590 
D  =  0.0 

DO  5  80  I  =  1 ,  a 

5  80  D  =  D  ( 7 (K,I )  /Yii  AX  (I)  )  *  *2 
?3  1  =  (D/£J«N)  **EN  Cl  *1.  3 

5  90  CONI  IN  J2 

IF  (PP2. L£.P3  3)  GO  TO  6  30 
IF  (P3  3.LT.  Pa  1)  GO  TO  660 

6  00  3  =  1.  J/.VUjil  (PS  1  ,  1.  E-4) 

N  2  ,tv  =  N0"1 
6  10  I  DOL'D  =  10 

I?  (  (KFLAO.  £o- 1) -AND.  (3.LI.  (1.  1)  J  )  GO  TO  700 
IF  (  NE*  *.  LE.  N  J)  GO  TO  630 
DO  o20  I  =  1 ,  N 

6  20  t  (  N2*  D*  1  ,  Ij  =  ZB30S  (I)  *  A  (  X)  /DFLGAT  (K) 
o  30  K  -  N  E  ’<i  ^  1 

IF  (KFuAG.  L*.  1)  GO  TO  670 
3ACJJ  =  a AO  JO  *3 
I BET  =  3 
GO  TO  7  i  0 

640  1?  ( NE«j.  £*.  .J^)  GO  TO  2  50 
N  0  =  N  2  »  2 
GO  TO  170 

650  IF  ( PE2.G1.  PL  1)  GO  TO  600 
NEa'G  = 

P.  =  1.  J/iUAXl  (?R  2  ,  1.  E-4) 

GO  TO  olO 

6  60  R  =  1.  J/AilA  X  1  (?ti  3  ,  1.  3-  4  ) 

N  E  W  v  =  N  ^  ♦  1 
GO  TO  610 
670  13  EX  =  2 

P  =  DM  IN  1  (5  ,  t(8  AX/D  A3  8'  ( H) ) 
ii  =  n«y 

LNErt  =  H 

IF  (  Ni* .  N  £  *  2)  GO  TO  6  30 

NO  -  NSNg 
GO  XO  1  7  J 
680  3  1  =  1.0 

DO  6  JO  J  =  2,  K 
21  =  31*6 
DO  6  1 0  I  =  1  ,  “! 

6  90  7  ( J  » 1 )  =  Y  (J,I)  *31 

IJOU3  =  K 

700  DO  7  10  I  =  J , N 

710  7.1  AX  (I  j  =  J.1i.(  1  (  YK  AX  (I)  ,  D  AoJ  (  i  ( 1  ,  I)  )  j 

0  OTA  3T  =  N  * 

F ZT  j  N N 

720  17  ( N 2 . E w -  1 1  GO  TO  7o0 

CAL-  DiFFJit  (T,Y,SAVE(M2,1)  ) 

3  =  H/riO-D 
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oo  7 .30  :  =  i,  :i 

i  t  1,1)  =  SAV  12  O  ,  I) 

GAVE(2,I)  =  :iGLD*5AYE  (NCI,  1) 

7  30  i  [2,1)  =  J«Vi(2  ,  I)  *3 

•'4  s  1 
K  F  LA  0  =  1 
00  I  0  1 7  0 
7  40  KFLAO  =  -1 
KNZrf  =  ii 
J  START  =  Nwi 
SETH  ON 

750  TUCGd  =  OlUXipABS  (HHIU/HOO)  ,3AC'J.1) 
RAC0H  =  iMIlil  (PAC'JH,  DA3S ( HH AX/HOLJ)  ) 
Si*  1-0 
DO  7  60  J  =  2,  X 
a  1  =  H  1  *a  AC 3 H 
00  7  o  0  1  =  1  ,  N 

760  7(0, 1)  *  SAVE  (J,I)  ■* S  1 

ii  =  iiOLO-»SACO;i 
60  7  70  I  =  1  ,  N’ 

7  70  1  (  1,1)  =  SAVE  (1,1) 

I  DC  0  0  =  i\ 

00  TO  (100,250,640), 10 ET1 
7  30  KFLAG  -  -4 
GO  TO  47o 

*j»t  0 


SJBBOUTIHE  DIFF'JS  (T,  77,077) 

IMPLICIT  -i  Z  AL*  3  ( A-  H,  j-  Z ) 

DICE  NS  UN  0  77  (23)  ,  77  (3,20)  ,  A  (20,20)  ,3  (20,23)  ,C  (20)  ,.(2(1)  ,»  A  (3)  , 
1  ^  A3  (20) 


0  *  *  *  *  *  *  : 


sgiuotinz  Called  37  cifsob  evaluates  the  derivatives  u? 

D  Z  ?  Eli  0  23  I  VARIABLES  STORED  1.1  77(1  ,1)  ,1=1,2*!.,  ..  N  3  STORES 
THE  DERIVATIVES  111  THE  A3HA7  077.  I  IS  THE  Ill  LEP  Ell  0  Ell  T 


/  A  i\  X  A  3  U  £t  a 


1  a  T  -  a  N  A  I*  A  J  X  1 

CO;iKwN/u1/COuF1#COEF5 

Cjr..10:i/L2/J0272 

w  3  /« kl  0 N /  u  3  /  .i  t  N  T 

C 2 .1 .1  ON  /L 4 /A  ,  o 

C  J . ifcl  J N /  u  7  /  i»  j*j  i  !l  2.  -  *  »l  ? 

ZL  =  2  *:i 
J  3  0  J.  =  1  /  1  I 

30  v.  j  »n  ^  N  J  Ci 


0  A  mw  j  0 a  r  0 


H  -  /A  w  Jwo  ^  N  ^  .i  j-  3  J  «?  .*  A  L  —  i  a  *  ) ) 


15^ 

ii(i)=YTli,  1) 

i  r  zs  n=  5  ou 

‘A  mU  Ll  *J  i,  ~J  X  »  i  ^  J  .C  1  /  H  r*  ,  |  i  l  Jf  ,  |4  *\  ,  ^  .1  U  ,  A,  'It  3  ) 

CT=-  J.  *CoZ7  1  *  (A  (1  ,  1)  *Lk  (1  )  +  A  I  1  ,  j  r )  *  1 .  ;  /  A  (1  ,  1}  -  A  (  1,  .,2)  /«  i  1  ,  1;  ‘U 
DO  *4  o  i  -  1  ^  d 

C  (I  j  =  3(1+1,  1)  *CT  +  5  (1+  1  ,NT)  +COZ72 
45  CONTINUE 

JTORi  IHE  D  E  F IV  A  T IV  Z  J  IN  AsRiii  D'{'{  — - 

30  1  20  J  =  1  , 

DVi  (J)  =  j  (0  +  1  ,  1 )  *.{.(  ( 1)  +  3  (3+  1  ,  NT)  *1  . 

3-Y  (J  +  N)  =0  iJ) 

33  110  K=  1  ,  :j 

DIi  l  J  )  =  3  1 1  ( J )  +3  ( J  +  1  ,  X  +  1 )  *  i'  i  (  1 ,  a  ) 

3  YY  (J+i<)  =  3  YY  (0+  N)  -3  -  *3  ( J  +  1  ,  1)  *C0E7  1  *A  ( 1  ,X+  1)  *YY  ( 1 ,  A)  /  A  (  1 ,  1) 

1  +  (3(J  +  1,  X+1  )-3  (0  +  1  ,  l)  *A  (1  ,X  +  1)/A  (1  ,  1)  )  *YY  (  1,.(+.<) 

110  CJ  NT IN  JO 

DYY  (J)  =302?  1  *C  0  E  ?  5  *  D  Y  Y  (J) 

Jii  (J  +  N)  =  J  u  (J+N)  *C02  F 5 
120  CON  i  IN  >1 L 
F  IT'J  SN 
EN  3 


j'J-JROJfiN  Z  PS3ZR  V  (I,  i,  ?  N  J ) 

IMPLICIT  SEAL*  3  ( A  -  H,  Q-  Z) 

DIMENSION  Y (3,20)  ,  P  W  ( 4  00 )  ,  A  ( 2  0 , 2  J )  ,3  (20,20) 


■-  *  a  : 


*  *  V  *  *  *  *  ■#  *  *  =  *  < 


:  **  *  *  *  Jt  * 


C  EJSRoUTINE  LAL-23  3Y  31: 003  STORES  DIE  PARTIAL  0 
C  LZ  NI7ATI /Z3  jF  THE  3IFFZ3ENIAL  Z2JA1ION3  PPOVIjSJ  In  C 
C  SUBROUTINE  Dlf’JON.  THE  PARTIAL  DERIVATIVES  (TOE  JACosiAJ)  C 
C  -ZPE  EVALUATED  IN  THE  LA5T  SECTION  IN  CHAPTER  5. 

C*************  =  *******»**»»****^=***=Ji***i!t«S*****#4i**»***!»*ii**«ji 


COMM  ON/ L 1/COE? 1 , CO  EF5 
CO.WON/LJ/N  ,  NT 
CC'IMoN/Ln/A  ,  3 


3o 

1  0 

i\  X  =  1 

,  J 

30 

4.C 

11  = 

1  , 

P .» 

l  II 

+  (2* 

r£  :i“ 

i)  * 

■l) 

=C  CL 

7  1  *Cu.j?5  * o  (H  +  1  ,NZ  +  1) 

L  •< 

(II 

*  (-* 

K  “ 

i)  * 

N) 

—  *  C  C 

Z  f  5  *  j  .  *  0  (11+  1,  1)  *COZr  1*\ 

1  1  r  Sr.  +  1 ,  / A  v  ■  .  "  / 

?  N 

(II 

+  <.*.! 

*  ;j  ♦ 

u* 

<:< 

-2)  * 

i)  =0. 

j  p  .» 

(  II 

+  2  *N 

*  ♦- 

(2* 

K 

)  -  i  3  ( 1 1  *•  1  ,.<•<+  1  ;  -  -  \I  I  +  ! , 

1 )  *  A  (  1 ,  a  .i  +  1 ,  /  A  ( 

10  Cu  .J  2  LN  J  Z 
.',ET  j  N 
ZN3 


n  (i 


! 
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susaoun.iE  DiCo:i?  (nn,a,'JL,  J2) 

DI.1ENS  ID  .'i  A  ^20  ,20)  ,UL(2  0,  2  0)  ,  SCALES  (2  0  )  ,  IPS  (20) 


:  **  «* «*  *  : 


********  *************.2 


SU330UTINE  CALLED  BY  MA7INV  PEfiFCdJS  A  FI8ST  STAGE  CF 
GAUSSIAN  EL  I M  I  NATION  SOUTINE. 


C 


1 

2 

J 

4 

5 


1  0 
1  1 
12 

1  3 

1  4 


1  5 


CO  Mil  ON  IPS 
J  2=  1 
N=NN 

DO  5  1=1,  N 
I  PS  (I )  =  I 

SO  NN  ail  =  0.0 

DO  2  J  =  1  ,N 

UL  (I,  J)  =  A  (I,  J) 

I?  (SO*  NS;i-a  ES  <U  L  (1,  J) )  )  1,2,2 

aOWN&il  =  A3  S  (UL(I,J)  ) 

CO  NT  I  ;i  U  E 

IF  (B^rt.ia.1)  3,4,3 

S  CALLS  (1)  =  1.0/P  OW  N  KM 
GO  TO  j 

CALL  SING  11) 

J  2-~  1 

SCALES  (I)  =  0.  0 
CONTINUE 
Nttl  =  N-1 
DU  17  K  =  1  ,  N  :i  1 
BIG  =  0.0 
DO  11  I  =  K ,  N 
I?  =  IPS  (I) 

SIuE  =  A3S  (QL(I? ,  X)  )  *  SCALES  (I?) 
IF  (SIDE-BIG)  11,  1  1,  10 

j  .  j  —  a  i. u 

1DXPIV  =  I 

.Oil  *  I N  0  ** 

IF  (BIG)  13,12,  13 
CALL  SING(2) 

J2  =  -  1 

GO  TO  17 

IF  (IDIP  I V  - K )  1  4,  15,  14 

J  =  IPS  (K) 

IPS  (K)  =  IPS  (IOXPIV) 

(I  JiCPIV)  =  J 
K  P  =  IPS  (K) 

PIVOT  =  JL(X?,K) 

X  Pi  =  K>  1 
DO  1o  I  *  X ? 1 , M 
IP  =  IPS  (I) 

E.1  =  -  JL  {IP  ,  K)  /?!  70  2 
J  I*  ( I  ,  a  )  *  -S3 
do  lo  j  =  x?i,:i 
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JL(I?(J)  =  UL(I?,J)  i-  i.1*'JL  (HP,  J) 

1b  CONTINUE 
17  *wQ  N T  xN  J  a 
X  ?  =  I  P  S  i  N ) 

:?  19,18,19 

16  CALL  SI  NO  12) 

J2=-  1 
19  PET’J  S N 
END 


SUoKOUTiaS  SOLVE  (JIN  , ?UL,3,X) 

IMPLICIT  o SAL*8  (A-H,w-Z) 

DI.iaNSi.OJJ  PUL  (20,20)  ,  3  (2  J)  ,\  (2  J)  ,  IPS  (20) 


'*«**#*****#*, 


C  SJJ3&OUTIJIE  CALLED  32  DIF3U3  PEPrGPHS  THE  SECOND  STAGE  C 
C  (THE  EACH  SUbSrirJTICtl  PROCESS)  OF  T.iS  GAUSSIAN  C 
C  ELI  Hit!  AI  ION  SOUTINE. 


i  v  m  *  +  *  m*  i 


‘  **  *c 


COMMON  IPS 
N  =  JIN 
N?1  =  N  + 1 
IP  *  IPS  i  1) 

X  (1)  *  B  IIP) 

DO  2  I  *  2 ,  !i 
I  P  =  i?s  U) 

I  ;i  1  =  i-i 

s  u  ;i  =  o.o 

DO  1  J  =  1,131 

1  sua  =  sj.i  ♦  pol  (ip,j)»x(j) 

2  X(I)  =  a  (IP) -SUJ1 
I?  =  IPS  (ji) 

X(N)  =  X(.«)  /POL  (I?  ,N) 

DO  4  I  SACK  =  2 ,  N 
I  =  NP 1-ISACE 
IP  =  IPS  (I) 

I?1  =  1* ] 

S  JM  a  0.  0 
DO  J  J  =  I ?1, N 

3  3UJ1  =  S0.1  +  PUL  (IP,  J)  *.<  (J) 

4  X  (I)  =  (aU)-SUH) /PUL  (IP,  I) 

PSI'U  oN 


SUoSJUIlNE  SINS  (I*HY) 


'*«***«  i 


*«*«***«i 


Sa33DJTI.it  CA -LED  JY  D  EC  C  HP  INDICATES  I  HE  ESSOsS  IJJ 


1 


u  <J 


15? 

PERFORMING  IJi  GAUSSIAN  ELI/.  I  NATION  aOJTINI. 


v. 

C 


i**************** 


‘******** ****** *******0 


1'  FOR  M  A I  ( 54  .i  J  H  ATR  IX  WITH  Z  E  HO  ROW  IN  DECOMPOSE. 

12  F  0  AM  AT  (  0  4  jj  3  S  I  NS  UL  A  H  SAIM  a  Id  DECOMPOSE.  ZErtO  DIVIDE  Id  _>CLVL. 

:  Id  IKPRU7.  MATRIX  ID  NEARLY  G  *.  Li  G  U  L  A  ? , 


13 

FOLilAr 

(54  HO 

NO  CONVERSE 

sour  = 

3 

so  ro 

(1*2, 

J)  , 

IWH  Y 

1 

WRITE 

(dour 

,11) 

SO  TO 

10 

2 

WRITE 

(doOI 

,12) 

vjO  t  o 

10 

3 

WRIT  E 

i^uU  r 

*13) 

10 

RETURN 

u  li  0 

GU3AOU 

TINE 

HAT  I 

dV  (PW  ,  N 

DiilE  NS 

ION  ? 

a  (400)  ,  A (20 

*  **  ****** 

**  **  * 

*********** 

3UDAJJ7I 

N  E  CALLED 

JY  DIF 

‘  *  *  *  C 


c 


GAUSSIAN  E  L  i  M  I  a  AT  I  0  d  Li  GUT  Id  E  OF  THE  H  ATP  IX  STORED  Id 
THE  AREA!  ?W. 

******** *8 i 


DO  1  I  =  1 ,  d 
do  1  J  =  1 ,  :i 

1  A  (I,  J)  =  II *  {J-  1  )*:;3) 
CALL  DEC  O  M  P  ( N  ,  A ,  U  L  *  J  1) 

R ITU RN 
END 


jUwRJU1I.»L  CU  5  EdT  (  CD  ,  A  ,  Y) 
IMPLICIT  RaAL*3  (A-K,  0-  Z) 

DI/.E  d5  10  d  A  (20,20)  ,Y  (8,  20)  ,  XS  (2) 


C  3 1J3ROU  TINE  CALCULATES  THE  CURRENT  DZNSITi  H'i  ScUAIIod 

C  ( 5 o)  . 

v.  *"* 

V*  W 

c  C3=caa«z;«;  jl*.oI;y,  a^p/c  ;<.«*<:  : 

C  A-  OliC  jZ 1 A r  lu  CO  £  IC  I  IN  I  *.  ATS  I.C  ,  A*  C 

u  * = c  <j  j  %*  w  n  j  it  *i  x  *  o  j  s  a  r  loLuula  Ti  points*  ^ 


*********** c 


n  r.  n  n  n  n  n 
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C  o  M  M  *yN  /Z  .j  /  .t  f  ,t  I 
CuMHON/..b/  XS 
common/Lo/co^e-. 

CD=-CJ&;4*  (  A  1 1  ,  1)  *X5  (1 )  *A  (  1  ,  .•  I)  ) 
Du  bo  J  =  1 .  :i 

CD=  Cii-ZOH?  4*A(l,J  +  1)*i[  '  1  ,  J) 

50  CONTINUE 
I ETU  Rli 
END 


S  U  3 N  j  U  T I  N  E  CO  N  C  7.  N  ( C 1  , C  2  ,  N  X  ,  N  D ,  D  ET  A  ,  E 0  0 1 , 0 1 F 1 ,  l j 
IMPLICIT  A  E  AL  *3  (A-  r.,  0-  1 ) 

DIMENSION  Cl  l  60)  ,  C2  (60)  , ROOT  (2.))  ,  OlE  1  (20)  ,  Y  (  1,^0)  ,  13  va,  ,  XT.  Nl  P  ( 2U ) 


C 

c 


S'JSf.OJ  TINE  EVALJATES  THE  CO  NC  ENT  u  A  "I  0  N  AT  EACH  DESIRED 
POINT  IN  X  DIRECTION  EY  AN  IN  I ZN PO LA T 10 N  ROUTINE  JOINS 
THE  KN  Oil  .i  CONCENTRATIONS  AT  THE  COLLOCATION  POINTS. 

~  —  —  .A  i  .I)  N  ~“ 

C  1  ,  C2=  CD  NC  E.i  1  LA  TIONS  u?  NITRATE  AND  n'YDSOXY  IONS  , 

RES  P  ECT  I  /  ELY  ,  IN  D  I  HE  NS  ION  L  E~>o  UNIT. 

D  E  T  A  =  X  I N  T  E  o  i  A  ,  DIMENSIONLESS,  0  <C  D  E  T  A  <  1. 


L 

C 


c 

i « ;• 


CO MM  ON /L 2 /CO EE 2 
C0MM0N/L3/N , N T 
COMMON/Lb/Xu 
DO  3  0  J=1 , N  X 
ETA  =DETA«  ( J  -  1) 

C  A «.  L  I N  YEP  ( .<  D ,  N  T  ,  ET  A ,  ROOT,  Die  1 ,  X  a  N 1 ) 
Cl  ( J)  -  X  i  L  TP  (  1)  *XS  (1  )  ♦XINTP  (NT)  -1. 

C2  (J)  =  XiNT?  (  1)  *X5  (2)  +XIN  TP  (NT)  *CO  ET  2 
DO  7  0  .•;=  1  , .« 

-  1  (J)  =C  1  (C)  +XINT?  (X*  1)  *  Y  (1  ,  K) 

7  0  Co  (J)  =C2  io)  «-X  I  NT?  (K+  1)  *1(1,  h+.i) 

30  Cu  N  T  IN  U  E 
R  E  TU  s.N 
END 


1  ET  A'j  M  A,  J  S  i  A .i  T ; 

2.1  PC  IC  IT  R  E  .i  w  *  S  ( A  -  H,  O-  E ) 


*>  '  J  -  .*  J  u  •  ^  u  1  >1  K.  e  .-'i  ^  .  J  .  Hy>  j  P  u  .  1  .  e  >  .  '1  1  ,  yw  i  A  .«  0  .  .1  .  y  4  >l.«i 

A  .» .i  .‘*e*  .  e*  -3  *  .  v. .!  A  I  L  ii.\  .i  y.  J  J  .  .  0  T  H  w  i  y  C:>  eAJ  .  y  .*  ~ 

JpI^Ny  P  0  I  .  ■  i. 
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*«****«#*  ******* 


******  *******  -1 **********************»********♦**£ 


GA7E:>=CS 

IF  (OH.  i  Z  .  0  .  u  0  1 .  A  N  D  .  2  S  .  L  T »  J  .  0  1  J ) 
IF  (-  J.  G  E .  J.  U  1  O.  A  N  3  .  Z  5.  LT-  0.  10  0) 
Ir  (u  5 «  G  Z  •  0  -  1  J  0 )  3  u  5=  0 •  05 

ZS  =  4.S*0C-. 

FACro»l=wZ/SAVZ3 

CO  OF  J  =  C0  if  j/f  ACTOR 

CjZF 4  =  00  if  4 /FACTO  3 

CO  ZF  3  =  00  Zf  5 /F  A  CTO  2/?  AC  TO', 

otaj  =  d:a  j/ i o j . 
itaj/.i*.  oj  i  *  j  ix’ j 

D I  AO  .1 A  =  2  .  Q  *  0  2  A  J 
J  5TA  &T  =  0 
fe  :u  a  u 
end 


0 ij i  —  0  .  0  005 
3iS=0 .005 


2032 30  TINE  EIPANO ( Y, ROOT, JiF 1 ,  F  ACT 03,  NO) 

ILLICIT  i>iAu*d  (  A  -  H,  0-  C ) 

3 1 H  E  .<  5  i  3  .<  t  (o  ,20)  ,  ROOT  (20)  ,  31 F  1  (20 )  ,  Yi  AT  ?  *  2  0  ,•  ,  0  Y  ( zC  )  ,  ii  ( 2 ) 


'*********  i 


S  0  H300  TI  .4  Z  S7  A- OATES  THE  CO  NC  E  0  T  S  ATI  0  00  AT  TiiE  CGLiCCAIlON 

POINTS  IN  li  ii  .(EA  CO  OF  01 X  AT  2  COE  TO  THE  CHANGE  OF  THE 

5? LINE  POINT.  C 

C 

■******«*****w«*»*»*******  *ii***»»*»4**y***<.*»*»****tf,; 


COHO  ON  /L  2/C  Gi  f  2 
C 0 .'Li UN/LJ/N  ,  NT 
C0.1H0N/L3/:n 
00  7  0  I  =  1  ,  N 
ST  (II  =  7  '  1  ,  I) 

7  0  0  Y  ( I  ♦  :i )  =  £  (  1  , 1  f  N ) 

DO  100  1=1,:) 

Ri=aOJr  (!♦  1)  *  F  ACT  CR 

IF(aE.OT.I.O)  GO  TO  50 

C.iL  L  INCR/  ( N  D ,  N  T  ,  BE,20<3T,  0 1 F  1  ,YINTP) 

7  (1  ,I/=tIL  4-/(1)  *XS(1)  <■  1 1 ..  T  P  ( .1 T )  *1.0 
7  ( 1  ,I*N)  =  TINT  ?(  1)  *X3  (2)  ♦  Cl  ..TP  (NT)  *CJ2F2 
Ou  ",  )  r.  =  1  ,  .< 

i  (  1  ,  .  j  =  ;  v  1  #  I )  f  Y  IN  T?  ( 1<*  1  )  7  ( ) 

a  0  Y  1 1 , 1  *■  N  )  =  Y  { 1 , 1  *■  N)  Y  I  N  T  ?  v  a  *  1 )  *  J  7  ( X  *  N  ) 

Go  . u  100 
9  0  Y  (1  ,  i  j  =  1 .  J 

7  t  1  , 1  *  ii )  =  C  0  i  F  2 
100  CCNCIN.Ji 
Ml  URN 

£.  .«  U 


n  o 
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F  J  NOTION  A U i l  ( XX  ,  K , J A3  ) 

I.iOLICIT  5  S Ao*6  ( A-  H,  *-Z) 

D I  IE  NS  13  N  ti  (/0,2  0 )  ,5  (20,23;  ,  X  X  ( 1 )  ,  2  A3  (20) 


C 

c 


«***************< 


t- 


C  A  J  X 1  IS  THE  NAJZ  OF  THE  FUNCTION  CALLED  31  ZSTST.1  IN  C 

C  THE  5U  OK JU 1 1 N  a  DIFFUN  TO  FU  So  IS  a  THE  VALUES  OF  EQUATIONS  C 

C  (52)  AND  (53)  .  C 


:****«*’ 


CG.M  JN/Ll/COt?  1  ,CO  EF5 
COMi1UN/Ld/COEF2 
COU.ION/LJ/N  ,  NT 
C 0 ilil 0 N / L 4 / A,  d 

C3dHJN/Ld/-uEF3,CO  N 1 ,CON2, SANA, BETA 

2  A3  (20)  =-3.  *CCEF1*  (A  (1  ,  1)  *XX  (1)  *A  (  1,  NT)  *  1 . 0)  /A  (  1  ,  1 )  -a  (  1  ,  :.T )  / A  ( 1  ,  1 ) 
1  *C  J  2f  2 

DO  5  00  J  =  1  ,  N 

2Aa  (20)  =  j  AN  (20)  -  3.  ♦COE?  1*  A  ( 1  ,  0  +  1 )  *2  A  3  (J)  /  A  ( 1  ,  1)  -  A  ( 1 ,  J+  1) 

1  *  < A3  (J  ♦  N  )  /  A  { 1  ,  1 ) 

5  00  CONTINUE 

IF  U AS  (20)  . LT. 3. )  CAP (20)  =0- 

AUX1  =A  ( 1,  1)  *  XX  (1)  *  A  (  1,  NT)  *1  .  ►  (CON2*  (QAi  (20)  **  3  A  HA)  -CON  1  * 

1  (XX  (1)  *  *  3  ET  A )  J/COEF3 

DO  (jOO  J=1,j 

600  AUX1  =  AUX1  *«t  (  1,J*  1)  *QA5(J) 

3ETU3N 

END 


FUNCTION  A  J X2 ( XX, K ,3  AS) 

IuPLIwIT  3  w  A i.. *  d  ( A  —  ri,  0"“-*) 

DIMENSION  A  (20,20)  ,5(20, 20)  ,I(d,20)  ,XX  (1)  ,  L  A  d  (  1 ) 


■**  *** «*«* . 


C 

C  A  0  X  2  IS  THE  NA.1Z  OF  THE  FUNCTION  CALLED  o'i  ZSYETJ  la 

C  THE  CAIN  FLuGBAH  TO  FU3NI5H  T  rfZ  VALUES  OF  Z<*UATIjNS 

C  (52)  AND  ( 5d)  - 


C 

C 


;**«*»**ae*«<*«a***’ 


CL  .IX  ON  /L  1/COS?  1  ,CO  EF5 
C  0  S.'l  0  A/ L2/C  UE  F2 
C  3.1.10N/L  J/N  ,  NT 
C  0  OH  O  N  /  L  N  /  A  ,  u 

Co  .1.1  ON  /L  3/  -  O  E?  J  ,  CO  N  1  ,CO  N2  ,  O  AC  A  ,  BETA 
CJ...1  JN/ui/i 

3a A  (  1)  =- J.  *C„EF  1*  (A  (  1  ,  1  )  *Xa  (1,  ♦  A  (  1  ,NT)  *1.0) /A  (  1  ,  1)  -  A  (  1  ,  NT)  /A  i  '  ,  1  ) 
1  *Co  E  ?2 

30  5  30  J  =  1  ,  N 
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3AS  (  1)  =EA  ji  i  1)  -  3-  *COE?  1*A  (  1  .  J  *  1 )  *  l  V1  ,  J)  /.i  ( 1  ,  1  i  -  A  v  1  1) 

1  *1  (  1  *  J*’  N)  /A  (1,7) 

500  cociruus 

ic  (OrtS  (1)  .  i.r.  0.  )  3.'Bp)=0. 

A  U \2  =  A  I  1  ,  1 )  *  A  a  ( 1 )  *  A  (  1 ,  :J  I)  *  1  .  ♦  (CC:i2  *  (  3  Aii  ( 1)  *  *  J/»  JA)  -CO  I.  1  * 

1  U-(  (1)  »*3ETA)  ) / CO  Z  F3 

DO  6  JO  J  =  1  ,  ,i 

600  AlU2  =  A  Ji2  *A  (  *'L  (1,  J) 

BiT’JSa 

E:i  o 


lOL 


JIUMItlXMXMOtMttXMOtM  **«***********««**************«. J 

C  c 

C  THIS  ?  ii  0  o  R  A  M  oU  L7  ^3  z.Q  U  A  T  10  N  3  {  4  d )  ,  (  4  ■)  )  ,  AND  (52)  i  r.  R  0  J*j  ii  w 

C  (55)  «xTa  I'tiJ  INITIAL  CONDITIONS  USING  SPLINE  COLLOCATION  C 

C  METHOD  FOR  THE  CASE  CF  THE  2LECTR0L  ITT  CONTA  IN  I  NO  CADMIUM  C 

C  IONS.  IN  THIS  CASE,  THE  TWO  SUBROUTINES ,  DIF FUN  AND  PEDLGV,  C 

C  CALLED  BY  THE  INTEGRATION  SUBROUTINE,  DIFSU3,  HAVE  TC  BE  C 

C  CHANGED.  MEANWHILE,  THE  MAIN  P20G3AM  AS  WELL  AS  THE  CTa.i  C 

C  SUBROUTINES  STAY  THE  SATIS.  THE  CONCENTRATION  PROFILE  Or  C 

C  THE  CADMIUM  ION  :1AT  BE  OBTAINED  BY  THE  ELECTRONS JT3ALITY  C 

C  CONDITIO.,.  THE  CURRENT  DENSITY  DATA  ARE  CALCULATED  IN  IHe,  C 

C  SAME  WAY  BY  EQUATION  (56).  THE  HETEROGENEOUS  REACTION  RATE  C 

C  CONSTANTS  OBTAINED  BY  LAST  P3CG2AM  ARE  USED  TO  IDENTIFY 


C  THE  HD  .10  SEN  £0  US  REACTION  RATE  CONSTANT  BY  FITTING  THE  C 

C  CALCULATED  CURRENT  DENSITY  DATA  WITH  THE  £.\?  S3 1  HEN  TAL  C 

C  DATA-  C 

C  0 

c  -  ADDITIONAL  NOMENCLATURE  -  C 


C  S  M  A  L LK  =  ri  D  M  0  G  E  N  E  0  U S  REACTION  RATS  CONSTANT,  1/SEC.  C 

C  D5IGK-0IMENS1ON  LESS  HOMOGENEOUS  REACTION  RATE  CONSTANT, 

C  DB  i  0  K  =  3  M<t  LLX  *  D  MA  X  **  2/5.  2b  £ -5  .  I 

C  D  KS  ?=  D I M  L  N  3  1 0  N  L  ESS  SOLUBILITY  PRODUCT  OF  CADMIUM  HYDROXIDE  ,  C 
C  DKSP=2.0a-23/C1 3**2.  C 


C**  : 


.*«»********=<  r 


IMPLICIT  REAL  *8  (A-H,  Q-Z) 

DIMENSION  A  (2  0,20)  ,3  (20  ,20)  ,C  1  (60)  ,C2  ( bO)  ,  C  (20)  ,  Y  (o  ,20)  , 

1  0  IF  1  (20)  ,DIF2  (2  0)  ,DIF3  (20)  ,  HOOT  (20)  ,  7  ECT  (20)  ,  E  n  T.  -  a 

2  Y  MAX  (20)  ,SAV  E  (24, 20)  ,PW  (R00)  ,  X3  (2)  ,  X.<  { 1)  ,  (  J)  ,  GAR 

EXTE  RN  AL  A  J  XL 

Co  MM  ON  /  L 1  /  ~  U  E  F 1 , CO  ZF 5 
CC  '..MON  /L2/Z  Ot  F  2 
C  'ML*,  .  j  ,  N . 

C  .  aM  O.t  /  L  B  /  A  ,  L 
COM.1UN/L5/XS 


COMMON /Lo/C  OS  ? B 
CO  MM  ON /L  7  /.i  GIG,  NEE  ,  Z? 
COMMON/LB/EOEF3 ,CO N1 ,CO  N2 , G  AM A , BETA 
COMMON/LR / » 

CUMMON/L 1 0/0  3  13  K , D  KS ? 


INPUT  DATA  — 


•i  aA  0  (5,310) 
R  EAD  (5,320) 
R  .A D  i  5 , 3  3  0) 
2  EAD  (5 , 3  4  u ) 


..  D  ,  1 ,  0  ,  N 1  .  ,u,  3.|  .  .j..  a  ,  ,  3  .4 

V-1  3,  CEB,  DM  AX,  DTIME,  N  i  ,  EG 
^  A M A  ,  J  ETA  ,  C  0N2  ,CO.,  1  ,  G M  A L L K 
N  3  x  vj  ,  N  Z.  ,  E  ? 


.  A L  C ’J  L A T  E  Go  HE  C 0  .* G  T  A N  T . 


:  j 


?  R  j  0  5  a  M  ~ 


COE F 1  =  (1.  j02D-5)  /5.  26  0-5 
COE  F2=C2o/C  1  3 


n  r.  (  i  (  i 


r  — 


163 

COSr J=  l  1.  0  0<.D-5)  *C1  S/DOA  X/'DS 

0.32  F4  =2.  *  i  1.  0  D2D-5)  «C  1  3«^o4J7./J1-iAX/c  J 

0  S 4  j  X=  3  44  ft  L  L  *  D .N  .A  \  *  w*  A.*.  /  5.  2  o  ■/“  3 
DKS?=(F.D-2  3j  /(CIS**  J.  ) 

DX=D.1Aa/(NI-1) 

3 In j-  (5.  iSi-5)  4 DTI  MZ/D  JAX/DuAa 
TAL'1  A.X  =  (5.  2oD-5)  *T  IX  Ai/D.1  A  2/D  .1 A  .( 

SDTA  U=DTAU 
d  ta  =  i  -  /  <  :i  :v  - 1 ) 
nt=n  * o + n  i 

:i  s«=  2*4-i 

-  2ZT  INITIAL  DIMENSION-ESS  ri.iS-J.  - 

7\U~  0. 


J^4  Ikl  U  4.  4  Z 


maximum  and  .u;;i:iu.i  dimensionlls s 

i  1  .X  ^  I .« C  3  S  ME  NT  TO  :  :  J  S  Z  D  IN  T  H  Z  I N  T EG  ?  .A  1 1 0  .« 
SUBROUTINE  D  IF  S  U3 - 


d:.\j.u  =  .  do  1  *dta  j 
D TAJ  /.A  -  2  .  0  *  DT  A'J 
*  a  II  Z  (  o  ,  J  b  J  ) 


-  £  y  *  L  u  A  T I  THE  COLLOCATION  POI.'I  IS  (aCOTS)  OF  I.iS 

JACOBI  POLYNOMIAL  OF  OBOES  N,  AS  WELL  AS  IAS 
Fla  ST  AND  SECOND  DERIVATIVES  Of  ThZ  PoLi'NuMIA- 
A I  r  ri  £  FOOTS. - 

-  .i  L  lm  J  C  j  a  1  ( 4!  D  1 4O  .4  J  f  4  <  1f  rt  m  g  j  41 1  i)  I  r  D  -i.  r  2  ,  Dir  i  f  j  •>  o  4 ) 

<<  .j  i  r  s  (& ,  j  5  0)  ( ?.  0  0 1  {  j  j  ,0=1  ,»»T) 


ALOJLATZ  TiiZ  DISCS  ZTIZ  ATI  O.'l  COEFFICIS.il  OAT! 


ID 

20 


10=1 

DO  <.  J  1=1/..  I 

C  A  -  L  0  :  0  F  A  ^  N  0  ,  N  ,  NO  ,  N  1  , 1  ,  I D  ,  0 1  ?  1  ,  01 F  2  *  D  if  .1 ,  SO OT  ,  V  -Cl ) 

DO  ID  J  =  1  , N  T 
A(I,J)  =  7-CT(J) 
wul.Tl  .i b  z. 

"™  ””  —  —  ^  j.Co  L(l*  i  4  4k  44  JiJ.44>4i4  4.  ..lJkl  DD^i4  4  IwiJk.l  4  4  il  4k'  4.  k'k  2  — 4.— - 


30 

*0 


DO  4=1.  .u 

Cfl  .*4  J  l'  k'  ^  -  *  -4  f  kl  f  4.4/  f  .4  I  f  4.  f  4  J  f  J  1  1'  lfJlk*^/j4  4  3  f  4. 

D^y  JO  u  —  I  /  ki  . 

D  (I,  J,  =  7ZCT  (J) 

j.,  i  n.i  j., 


1.1PJ.  . V  A  L  .J  a  S  Jr  AF-jjkik^.T^  4  j  .j  Z  JSc. 

1  N  4  _  j  3  \T  1 0  N  S  J  j  30  ■)  I  I.N  C  J...  jJj  —  —  — 


1 
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HF=1 

0  S  T  A  aT  =  0 
ISTD ?=0 
HAXD£3=7 
DO  47  I  =  1  ,  it  Eel 
47  Y.'lAX  (I)  =1.  J 

1  =  0 
n 

C  -  INPUT  THE  INITIAL  CONDITIONS  - 

c 

DO  50  0=1,11 

7  11  ,0) =1- 

50  7(1  ,J+N)  =00072 

XS  (1  )  =  1. 

XS  (2) =  C0Er  2 
WRIT E  (6 , J90  ) 

60  TIS£=TAU*DRAa*DHAX/(5.  2  6D-5) 

v-  ———————  wALk-ULATE  AMD  ?8IMT  THE  CURu— NT  DENSIT  L  AT  To..'!  — 

0  PRINT  THE  CONCENTRATIONS  0?  NITRATE  AND  HYDROXY  IONS 

C  AT  Tii E  NT  CO LLOCATI CM  POINTS - 


call  c u s  e n x  (c d ,  a  ,  y ) 

"RITE  (6  ,  ioU)  TINZ,CD,XS(1)  ,  XS  (2)  ,  X  F  .u  A  G ,  I X  S3  N  ,  AS 
NRI  IE  (6 , 70  Oj  ( Y  { 1  ,L)  ,  L=1  ,  N2R) 

C  -  EVALUATE  AND  PRINT  THE  CONCENTRATIONS  0? 

C  NITRATE  AND  HYDROXY  IONS  A  l  EACH  POINT 

C  IN  X  DIRECTION  AT  THIS  TIME - 

C 

CALL  EON JEN  (C 1,C2, NX,ND,D ETA, ROOT, DIE  1  ,  Y) 

WRITS  (6,200)  (Cl  (J)  ,J=  1,  NX) 

WRITE  (o,250)  (C2  (J)  ,0=1, NX) 

C 

C  - CHECK  IF  THE  INTEGRATION  TIHE  AT  THIS  HU. 12  NT  REACHES 

C  THE  HAXillUa  TINE  TO  3E  INTEGRATED - 

C 

IF  (IAU.  GE.  TAUMAX)  GO  TO  190 
35  1  =  1+  1 

IF  (I.  EQ.  1 )  GO  TO  9  0 
IF  (ISTDP.  E«i-  1)  GO  TO  90 
C 


n 

c 

c 


TESr  IF  THE  CONCENT  RATI  ON  0?  NITRATE  ION  AT  THE  NTH 
E  0  u  LCCA  TI 0  N  POINT, IE,  T.iS  LAST  ONE  3EFORZ  THE  SPLINE 
P C  *  a  I  IS  *  IT  HI  N  J-  OOC  1  Dl.i  ZNS  I  JN  LESS  C  D  NCEN  IR  . I  w  i< 
UNI  I  OF  THE  30UNDARY  CONOITION.  - 


C 


C 

C 


I?  ( (  1.  -7  ( 1  ,  5)  )  .  LT.  0.  0C0  1)  GO  IG  )0 


IF  THE  TEST  FAILS,  THE  SPLINE  POINT  L3  I  NCI  LAS  E  J  ,  IE, 
IJVi„  FURTHER  INTO  THE  SOLUTION.  - 


V. 

C 


165 

lAuL  w  ri  A  N  vj  i  [  uS  f  Du  Sf  L.C  x. 3  t  ■-  JuC  A  |  D  /  £A'.I03,DTAU,D.AOHl,DxA  J  F.  t\  , 

J5TA3T) 


■~— —  nVALJATE  T  iiE  CO  MC  £..«  13  AT  ION  J  J?  N  IT  ?  A  7  E  ANO  HYOSOXY 

IONS  AT  THE  COLLOCATION  POINTS  AT  THE  N2X  CCGtiOI.UrX 
Oil  2  TO  THE  INC  3  EASE  0?  THE  SPLINE  POINT. - 


CALL  EXPAN  0  ( Y  , SOOT, DI F 1 , FACTOR , NO) 


- INI  EG  3  AT  E  EQUATIONS  ('48)  AND  (4j)  USING  THE  C  J  NC  E  N  13  AT  ION: 

AT  In 2  COLLOCATION  POINTS  AT  THIS  TI.iE  TO  OuTAaN  I  .:E 
CONCENTS  AT  IONS  AT  THE  SUCCESSIVE  TINE  - 


90  CALL  DIFSU3  ( N  2  It  ,  TA  U,  Y, S  AYE ,  OTA  J ,0 TAUH I  ,  OTA  UNA  ,  EPS,  NT,  l  A  A  X  , 

1  £.3303,  KFLAG,JSTA3r,HAX0  S3  ,  ?  A ) 

IF( XFLAG. LT. 0)  GO  TO  190 
DTA UMI=. 00 1 *0TA0 
OTA  UM  A  =  2. 0  *D  TArJ 
IF  (I.  GE.  2)  Go  TO  93 

IF  (Y  (1,1)  .0  2.  1..0R.Y  (1  ,2)  -  IT.  1 . . 03 . Y  (  1,  3)  .  GT .  1 .  .  O.i .  i  (1  ,4)  .  ST.  1  . 

1  .03.  DAoS  (  Y  (  1 , 5)  -  1.  )  .  GT.  0.  0004)  GO  TO  90 

GO  TO  05 
93  .11  =  1/10 
22=ai*IC 

IF  (I.EQ..12.03.TAU.GE.TA  0.1  A  X )  GO  TO  95 
GO  TO  d5 

- CALCULATE  THE  S 'J?  FACE  CONCENTRATIONS - 

95  OLDXS2  =  XS (2 ) 

XX(1)=Y(1,  1) 

ITL3  N=500 

CALL  ESYSTS  (aUX2 , E P, NS IG, NE Z, XX , IT E3 N , - A , 3 A 3 , I ER) 

XS(1)  =XX  (1) 

XS(2)  =-i.  *COE?1*  (A  ( 1,  1 )  *XX  (  1)  4  A  (  1,  NT)  *  1 . 0 )  /A  i  1  ,  1 )  -  A  (  1  ,  N  T )  /A  (1  , 1 ) 
1  *C0EF2 

00  150  J  =  1  ,  N 

XS  (2)  =XS  ( A)  -  3.  *CCEF1  *A(  1 ,  J*  1)  *Y  (  1  ,  J)  /  A  ( 1  ,  1)  -  A  ( 1  ,  J  ♦  1) 

1  *i  (1,  J«-N)/A(1,  1) 

150  CONTINUE 

IF  (XS  (2)  .  LI. OLD X 52)  ISTO?=1 
GO  TO  oO 

190  X3IT  E  (a  ,  2  40  )  D(,3  OTAU,  NX,  OTAJ,  DTI  HE,  OH  AX,  AL,SZ,  ES,  EPS,  .<  F  EG  j, 

1  E? , N  SIS 

.  SET  A,  GA  M  A  ,  CO  N  1  ,  C  J  N  2 ,  S  .1 A  E  L  \ 

- F  03  .UTS  FOR  INPUT  A  NO  UJYPJT  ST  A  2  El  EN  TS - 

230  FO  AN  AI  (/,2  5  X,  *  NO  3-  ',  IX,  :.D  11. 4/3  OX,  90  1  1  .  4/3  OX,  sJ  1  1 . 4/ j  OX,  30  11 
1  /  iO  a  ,  90  1  1  .  4/3  OX  ,  9D  11  .4/3  OX,  301  1 . 4/3  0  X  ,  9  s  i  1.  4/ j.}  ,v  ,  3D  1 1  .  4 


.i 


* 
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2 

2  40 

1 

2 

3 

4 

5 

2  45 

1 

2 

3 

2  50 

1 

2 

3  10~ 
320 
330 

3  40 
3  50 
360 
3  70 
300 
3  00 
7  00 


/JO  2  ,  00  11 .  4/302,  OD  1  1 . 4/3  02,  30  11 .4/302, 9  0  11.  4/ JO  2,  Ju  11. 4) 
FOBMAi  (//,5X, 'X  INTE3VAL='  ,  012-  4,  32*,  '  INITIAL  IAU  STEP  SUE=' 

,0  12.4/52,  '  SO  Or'  POINTS  10  2  013=  ',13,2  92, 

•LAST  TAT  STEP  SUE=  '  ,01  2.  4/ 52, 'PEAL  TiaE  INTSB7 AL=  ' 
,010.4,242,  'DIFJSION  LENGTH3  •  ,010.  4//  52,  '  A  L  3  ',.-7.2,302 

'  0  li  =  •  ,F7.  2,252  ,' ZJ=  '  ,  F7 . 4//5  X  ,  '  S?3=  ',010.2,27a, 

'  K  FL  AO=  '  ,  I3//52,  '  2?  =  '  ,  J  1 0 . 2 , 27  X , '  :i  3 10=  ',13) 

POa.i  AT  (//,  52,  'REACTION  03D2H  OF  SC3-  ION  -'.cl.  2,25a, 

'  d  2  ACT  10  N  QSDE3  OF  OH-  UN  =  '  ,  F  7 . 2//S  a  ,  '  1ST  CONSTANT  =' 

,  0  12.4,252,  '2ND  CONSTANT  =',012.4 
//52  , ' 3 2 ACTION  RATS  CONSTANT3  ',012.4) 

FOaa  AT (/ , 25  2 , 'OH-'  ,2X,  901  1.  4/ 30a, 901  1.4/302,001  1.4/302,901  1.4 
/JO  a, OD  1  1.  4/3  02  , SO  1 1.4/3  0  2,  901  1.4/3  0  2,9  01 1.  4/  jO  2 , 9  0  1 1  .  4 
/3  02  ,  901  1.  4/3  OX,  90  1  1  .  4/3  02,  9  01  1.4/33  2.  93  1  1.4/30  2,  ih  11.4) 
FOSHAT  {413, 2/5.  1,  2010.  2,1  5,F3.4) 


FOR:1AT(401  1.4,I7,F7.4) 

FORMAT  kzF7.  2,3012.  4) 

FORMAT  (215, 012.4) 

FORMAT (152, 1  OF  10. 6) 

FORM  AT  (//,  j  12. 4, 1  2  ,D  12.  4,  22 ,2030.  10, 215,  F20.  7) 
FORM  AT (/,  2  2, qF20.  6/2  X, 6:20. o) 

FORM  AT  (2  2,  '  COLLOCATION  POINTS:  ') 

FOR  a  AT  (//,  j  A.  'TI.7S  '  ,  12  A  ,'C0.Ha2NT'  ,45  X,  '  ION  CON 
FCSM  AT  (302,  1 0  F 1 0  .  5/30X,  10F10.  5) 

STOP 


ENTS All ON  '  ) 


END 


5  U330U  TINS  01 FFON  (T,  YY  ,  DY  Y) 

IMPLICIT  RSAi.*>3  (A-H.Q-E) 

01  a£  NS  ION  0  YY  (20)  ,  YY  (B,  20)  ,  A  (20, 2  0)  ,3(20, 20),  C  (20)  ,  22(1)  ,  «A  (3)  , 

1  ^ An  (20) 


C  SUBROUTINE  TALL  2D  3Y  0IF303  EVALUATES  THE  DERIVATIVES  OF 
C  DEPENDENT  VARIABLES  ST  03  EO  IN  YY  (1,1)  ,1=1,2*11,  AMO  STORES 
C  THE  DERIVATIVES  IN  THE  ASSAY  DU.  T  IS  THE  INDEPENDENT 
C  VARIABLE. 


EXTERN  «*L  A ;J  21 
COaaON/L 1 /C OEF 1 , CO  EF5 
Cuaa JN/L2/COEF2 
coaaoN/L J/j , n  t 
c  j  aa  ON  /  L 4  /  a  ,  3 
COa.I  ON  /u  7/.«  S.  j  ,  N  r.  E  ,  o? 
CoaaoN/LIO/ObIGK, OKS  ? 
I 1=2 *H 

JO  3  0  1=  1  , 1 1 
*  AS  (I)  =  Y  Y  (1,1) 

30  CONI IN JE 


c 
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CALCULATE  THE  7ALJZS 


J  U  .1  4  IlC  w 


V  A=0  > 


04 


. 4l  . 


4‘Yl1)=YY(1,  1) 

I T  E  B  4  =  3  u  0 

v.  t\  L  u  uo  /  j  i  ik  ( A  U  a  1  ,  --  ?  /  a  S  a.  kj  ,  aA,^Ab,^lb) 

CT=~3.  *CGEF  1*  (A  (1  ,  1)  *XX  (1)  *A  (1  ,  .i  r )  *1.  )/A  0 , 1)  -  A  l  1  #  NT)  /A  i  1  ,  1)  *CCEF2 
DO  4 3  1=1,4 

C(I)  =  u  ( I  *  1  ,  1)  *CT  ■*•  B  {  !•*•  1  ,41)  *C0£F2 
45  CONTINUE 

- SUEZ  THE  DERIVATIVES  I A  A  A  4  A  Y  DYY - 

DO  1 20  J  =  1  ,  4 

DYY  (J)  =o  ( J  *  1,1)  *  X  i  ( 1)  *3(  J+  1  ,  .«Tj  *  1  . 

DYY  (J*4)  =0  (J) 

DO  110  K  =  1  ,  ;; 

D  YY  (J)  =  DYY  (J)  +3  (J*1,  X*  1)  *YY  (  1 ,  K) 

110  D  YY  ( J  ♦  N  >  =0  YY  (J*  4)  -3  .  *3  { J  *•  1  ,  1)  *CG£F  1  *A  (  1  ,K*  1)  *  Y  Y  (  1 ,  a  )  /  A  l  1 ,  1) 

1  *  ( 3  ( J  *  1 ,  K  *  1 )  -  3  { J  ♦  1 ,  1 )  *  A  ( 1  ,  K  ♦  1 )  /  A  i  1  ,  1  j  )  *  Y  Y  (  1  ,  «v  ♦  4  j 

DYY  (J)  =G0£F1  *C CE F  5*  DY  Y  (J  ) 

DYY  ( J  «■  4 )  =0  YY  ( J  >  N  )  *C  OE  F  5 

SO?  £3=2.  *  0  S I G  K  *  (0 .5*  (YY  { 1,  J)  *Y  Y  ( 1 ,  J*-  4)  )  -DXSP/Y  Y  ( 1  ,  JO)  /Y  Y  .1  ,0  ♦  1.)  ) 
IF  ( S  li  ?  E  A  .  3  1.  0.)  DYY  (J  «■  4)  =  D  Y  Y  (J»4)  -SUPER. 

120  CJSTIUJc 
it  l  T  0  b  4 


S 036 OUT 14  L  2ZDE3 V  (T, 7,  ?4, 43) 

IMPLICIT  B  E  AL*8  ( A  -  H,  Q-Z  ) 

DIME4S1J4  Y (3, 20)  ,?W  (4  00)  , A  (2  0,20)  ,3  (2  0,20) 


S  0  330 J  TI 4  £  CALLED  3Y  DIPSOS  STORES  THE  PARTIAL 
DERIVATIVES  DP  THE  0  IF  FEW  EMI  A  L  STATIONS  PbOVIDED  I., 

S  'J2HOU  TI4  £  THE  PARTIAL  DERIVATIVES  (THE  JAC0SLA4) 

4ZF.E  EVALUATED  14  THE  LAST  SECT 104  14  CSnPTE?  5. 


C 

C  *  * 1 


U 

c 

c 


k  «  *  *  *  »*  i 


!  **  *  *  *  *  =  k 


COMM  OU/L 1 /C02F 1 , CO  EF5 
COMM  04 /L J/H  ,  N T 

Ckjktu  Ukl  /  . ■*  /  kk  ,  3 

C  CM .1  04  /l  1  0/0  a  I IX  ,  D  XS  ? 

Du  Is  rvK=  1  ,  A 
DO  20  II=1,j 

J*  (II  +  (2*  ab-  2)  *  4)  =CCEF  1*CC£F5*d  (II  «■  i  ,  aK  ♦  1 ) 

P  »  ( II  ♦  l  2  *  E  a -  1 )  *  4  )  =-CC£?5*J  .  *  3  ( 1 1  ♦  1  ,  1 )  * C  0  E  r  1  *  A  (  1  ,  X;.  *  1 )  /  A  k  '  ,  1  ; 

IF  l  II.  Z*.  K£)  GO  TC  30 

Jkj  k  D  4y 

30  Sl?=2.*JjT0K*(G.5*(Y(1,II)  ♦  Y  (  1  ,  1 1  ♦  4  )  )  -  D  .\ «,  ?/  Y  v  1  ,  1 1  ♦  4  j  /  i  k  1  ,  j 

IF(JU?.UT.O.)  P  *  ( II  ♦  (  2  *<K-  1)  )  =  2  *  ( 1 1  *■  (  2  *  £  K  -  1 )  *  :. )  -  OS  I  G  v. 

PA  (  II*2»a*4*  ( 2  *  K  a  -  2  )  *4)  =0. 


40 


1 
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?*  ( II  0  *  H  *  »<  ♦  (2*  XK  —  1)*M)-(.JiII*1,i\.;*'1)-b(II*1/1)  *A(1,Ki\+1)/Al1f1)) 

♦C0SF5 

Ir(II.i^.KjC)  JO  TC  u  0 
JO  TO  20 

60  IF(  jjp.  ji.  o. )  ?w  (n*2*M*.Sf  {2*K\-  ij  *:i)  =?*  u*.\K- 1)  *  jj  -2.  * 

1  0313  X  *  (.  6  *-2.  *0ao  ?/T  (  1,  II  +:• )  / 

2  i  ( i ,  1 1  +  :i )  /  y  ( i ,  1 1  *•  j  )  j 

20  CUSTI.iJJ 
10  CO. IT  IN  US 
l\  2  T  J  d  N 
END 


***************, , 


»******♦*«(; 


THE  CTrlEii  J  J  3  30  'IT  I tf  E3  CAL1.20  J'{  THU  PSuJTAH  STAY  T.iE  SA.1E  C 

as  those  ::i  the  pisst  peoosah.  c 


****** *«*************r 


* 


